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1. 


INTRODUCTION 


During  re-entry  into  the  atmosphere  of  a  hypersonic  ballistic  missile  the  char¬ 
acter  of  a  radar  return  changes  in  several  respects  from  that  of  the  return  observed 
when  the  vehicle  is  in  the  exoatmosphere.  Many  of  these  differences  (in  particular,  the 
often  rapid  and  random  variations  in  amplitude  and  phase)  are  believed  to  be  due  to  the 
presence  of  a  turbulent  wake  acting  as  a  statistical  scattering  medium.  It  is  generally 
assumed  that  there  exists  a  random  distribution  of  turbulent  eddies  which  are  large 
enough  in  number,  within  one  resolution  volume  of  the  viewing  radar,  to  permit  appli¬ 
cation  of  the  Central  Limit  Theorem.  In  other  words,  the  resultant  scattered  field  is 
due  to  the  superposition  of  a  large  number  of  components  with  independent  phases. 

This  then  implies  that  the  resultant  quadrature,  or  the  Cartesian  components  of  the 
scattered  field,  form  a  pair  of  independent  Gaussian  stochastic  process.  Thus,  the 
amplitudes  of  the  received  signal  are  Rayleigh  distributed,  i.  e.  the  received  power  has 
an  exponential  distribution. 

Often  the  resultant  signal  at  the  radar  is  not  that  due  to  the  turbulent  wake  alone, 
but  is  a  superposition  of  this  component  with  the  coherent  return  from  a  hard  body 
scatterer.  Such  a  superposition  may  be  observed  in  the  range  cell  containing  the  hard 
body  which,  because  of  the  size  of  the  resolution  volume,  also  contains  some  of  the 
turbulent  wake.  In  the  case  of  an  extended  return,  the  same  situation  may  obtain  even 
when  the  observations  are  not  made  on  the  range  cell  containing  the  body.  This  is  due 
to  range,  and  perhaps  angle,  sidelobes  of  the  body  return. 

In  this  report  the  statistics  of  the  signal  received  when  the  hard  body  and  wake 
returns  are  combined  are  investigated.  Four  computationally  feasible  methods  for 
estimating  the  component  powers, arising  separately  from  the  turbulent  wake  and  the 
hard  body,  are  developed  and  compared.  The  first  two  methods  utilize  amplitude 
data  only,  and  thus  are  suitable  for  general  use  with  the  data  obtained  by  most  radars. 
The  third  and  fourth  methods  utilize  phase  returns  as  well  as  amplitudes,  and  thus 
presuppose  a  coherent  radar.  An  intermediate  step  in  the  third  method  also  yields 
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an  estimate  of  Doppler  frequency. 

Specifically,  then,  the  raw  data  for  a  coherent  radar  are  pulse-to-pulse  am¬ 
plitudes  and  phases  within  a  series  of  range  cells.  In  accordance  with  the  preceding  dis¬ 
cussion,  the  model  hypothesized  for  the  amplitude  r  and  phase  p  recorded  on  the 
th  t  t 

t  pulse  is,  when  written  in  Cartesian  coordinates, 

r  cos  p  =  x  =  pcosw  +  <rw  . 
t  Kt  t  t  It 

(t  =  1,  2, . . . ,  n).  (i.  1) 

rsin  pt  =  y  =  p  sin  ^  +  crw2t 

(The  notation  does  not  indicate  the  fixed  range  cell.)  {w  }  and  {w  }  denote  inde¬ 
pendent  Gaussian  processes,  each  of  which  is  presumed  white  with  mean  0  and  vari¬ 
ance  1.  For  the  present,  a  linear  phase  model  is  postulated: 

cp ^  =  k>t  +  (p.  (1.  2) 

The  immediate  problem  is  joint  estimation  of  the  wake  and  hard  body  cross  section 
parameters 
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1 


(1.3) 


without  knowledge  of  the  phase  parameters.  We  are  chiefly  interested  in  rapidly 
computable  estimates,  i.  e.  those  having  storage  requirements  which  do  not  depend  on 
the  sample  size  n. 

All  evaluations  will  be  made  on  the  basis  of  results  of  large  sample  distribution 
theory.  In  practice,  n  cannot  be  taken  too  large,  because  the  assumption  of  a  con¬ 
stant  hard  body  amplitude  would  most  certainly  be  unrealistic.  Nonetheless, 
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asymptotic  comparisons  have  qualitative  comparative  value  when  considering  fixed 
sample  size  applications.  In  addition,  the  labor  required  to  obtain  higher  order 
approximations  to  mean  square  errors  is  not  worth  the  effort  when  the  validity  of  the 
model  is  still  open  to  question.  In  particular,  there  is  empirical  indication  that  the 
regressions  might  be  superpositions  of  sinusoids,  or  that  the  wake  processes  are  not 

t 

white.  These  questions  are  deferred  until  a  later  report. 

Since  there  may  be  an  additional  acceleration  term  in  (1.  2),  it  is  customary 


(for  incoherent  radars,  necessary)  to  focus  attention  on  the  amplitudes  r  =  v  x  +  y  , 

which  are  distributed  the  same  for  every  t  and  independently  of  the  phase  model.  The 

common  variate,  say  r,  is  said  to  be  Rice  distributed  when  p  f  0  (see  p.  238  in 

Wax  [7]),  and  Rayleigh  distributed  in  the  special  case  p  =  0.  In  statistical  terminology, 

r/or  has  a  noncentral  chi  distribution  with  2  degrees  of  freedom  (written  y'})  depending 

only  on  the  parameter  p.  When  p  is  0,  the  adjective  "noncentral"  is  dropped  and  the 

prime  deleted  from  the  chi.  This  distribution  obtains  for  any  two  mean  value 

2 

functions  whose  sum  of  squares  is  p~  (see  e.g.  Kendall  and  Stuart  [4  ],  p.  227).  There 
is  a  certain  notational  ease  inherent  in  the  statistical  terminology  which  makes  it 
worth  while  using  in  the  course  of  analysis. 

2.  SUMMARY  OF  THE  METHODS 

Throughout  we  will  let 


n 


(2.  1) 


t=l 


respectively  denote  the  sample  and  population  amplitude  moments  of  integer  order  k. 
The  dependence  of  the  former  on  n  is  not  shown  and  must  be  remembered. 

Method  I  is  based  on  the  fact  that 


t  The  results  which  we  derive  in  Sec.  7,  and  specialize  to  the  analysis  of  the  third 
and  fourth  methods,  apply  to  more  general  models. 
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e 


A 

2 


=  F(A) 


a 

1 


i 


(1+A)1 


where  1^  is  the  modified  Bessel  function  of  the  first  kind  of  order  vy  and 


a  =  o2/e 


1 


(2.2) 


is  the  signal -to -noise  ratio  (or,  more  appropriately,  the  hard  body  -  to  -  wake  ratio). 
The  function  F  increases  monotonically  from  F(0)  =slir/2  to  F(°°)  =  1.  The  ratio 
(2.  2)  can  therefore  be  estimated  by 


A*  =  < 


F  V/ap 


if  a  J \Itt/2 
1  2 

otherwise  . 


(2.3) 


Then,  since 


&2  -  expected  received  power  =  9^  +  0,^  or  0^(1  +  A), 


the  quantities 


I. 

e*  -  _*2  . 

1  1  +  A* 

92=v0; 

(2.4) 


are  taken  as  estimates  of  the  parameters  (1.3).  Figure  1  is  included  to  provide  a 

i  -f 

means  for  roughly  estimating  A  from  the  observed  ratio  a^/a^j  • 


t  The  estimates  (2.  4)  have  been  considered  by  D.  Batman  of  MIT  Lincoln  Laboratory. 
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We  point  out  here  that  the  (or  any)  Maximum  Likelihood  estimates,  based  on 


amplitude  data,  satisfy 


Wa2 


We  will  invariably  use  this  as  one  of  our  estimator  equations. 

2 

Method  II  is  based  on  the  fact  that  the  moments  of  r  /0  are  just  Laguerre 
polynomials  in  -A.  As  a  consequence 


6 


2 


Using  the  first  two  sample  power  moments,  then,  we  immediately  have 


II. 

1 

\  (2a^  -a4)“  if  2a“  ^  a4 

0*  =  a0  -  0* 

0*  = 

12  2 

2 

0  otherwise 

(2.5) 


as  consistent  estimates  of  (0^,0^).  It  should  be  noted  that  both  pairs  (#*,<?*)  are 
necessarily  nonnegative. 

For  either  method  of  estimation,  the  distribution  of 

^n(e*-ei)  ^n(e*-e2) 

tends  to  a  zero  mean  bivariate  normal  distribution  as  n  —  00 .  In  each  case,  the 
covariance  matrix  is  of  the  form 

2=^Q(A)  (A  >  0). 

The  joint  asymptotic  efficiency  of  Method  II  relative  to  Method  I  is  customarily 
measured  by  the  ratio  of  determinants  (called  generalized  variances) 
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(2.6) 


where  the  meaning  of  affixes  is  clear.  In  spite  of  the  fact  that  Method  II  is  computa¬ 
tionally  simpler  than  Method  I,  since  the  latter  requires  table  look-up  and  interpola- 
_  \ 

tion  for  F  ,  the  function  (2.  6)  exceeds  unity  throughout  some  neighborhood  of  A  =  0 
(below  some  point  just  short  of  .4).  Both  methods  become  asymptotically  efficient  in 
the  limit  of  indefinitely  large  X,  and  Method  II  also  does  as  X  —  0.  This  is  not  too 
interesting  since  the  generalized  variance  of  the  ML  estimates  tend  to  infinity  quad- 
ratically  as  X  —  0  and  linearly  as  X 

The  poor  statistical  behavior  of  (2.  4)  and  (2.  5)  for  small  values  of  X  can  be 
eliminated  when  we  use  the  phase  data,  and  assume  the  model  (1.  2)  holds  with  some 
arbitrary  frequency  and  initial  phase. 

The  time  series  in  (1.  1)  are  clearly  nonstationary.  They  are,  however, 
"asymptotically  stationary.  "  The  reason  for  this  terminology  is  that 

n- | h |  n- | h | 

A,  =  -  /  xx  i,i  B  =  —  /  y  y  ■  i  (lh|  ^  n-l) 

h  n  L  t  t+| h |  h  n  L  tJt+|h| 

t=l  t=l 

converge  as  n  —  °°,  both  with  probability  one  and  in  mean  square,  to  certain  limits. 

It  is  true  for  their  sum  that 

S,  =  +  Bl  —  0.  <5U  _  +  0  cos  h«  =  y  (-°°  <  w  <  <»,  X  0)  (2.  7) 

h  h  h  1  h,  0  2  h 

for  every  fixed  h,  where  6  is  Kronecker’s  delta.  The  Fourier  transform  of  this 
sequence  (h  =  0,±l,±2, . . . )  is  constant  at  all  angular  frequencies  save  at  to,  where  it 
has  a  jump  of  magnitude  9.2  over  and  above  the  level  9^.  This  fact  suggests  basing 
estimation  schemes  on  empirical  (windowed)  power  spectra.  Such,  however,  entail 
calculations  of  a  number  of  covariances  which  in  principle  increases  without  bound 
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with  sample  size.  For  the  present,  we  will  not  be  concerned  with  spectral  techniques. 

On  the  other  hand,  we  can  still  consistently  estimate  (0 ^  0 ,;)  by  various  simple 
combinations  of  the  first  few  covariances.  These  methods  can’t  attain  full  efficiency, 
but  there  is  an  improvement  in  variances  over  Methods  I  and  II  which  increases 
without  bound  as  A  gets  smaller.  Method  III  is  one  such  procedure  which  first 
estimates 

£  =  cos  cu  (2.  8) 


by 


s2  +  Js;  +  8S2 


4S 


(2.9) 


This  is  derived  from  the  fact  that{xj  and  {yj  obey  the  same  second  order  stochastic 

difference  equation  driven  by  independent  second  order  moving  averages  with  common 

variance  i9,(l  +  2£  ).  (2.  9)  is  obtained  as  the  (statistically  equivalent)  unambiguous 

1 

point  minimizing 


1 


1  +  2£ 


2 


y  [<xt  -  2|  xt_4  +  xtl)2  +  (yt  -  2?  yt.,+  yt)2l  • 
t=3 


(2.  10) 


This  particular  criterion  ignores  the  first  two  correlations  of  the  driving  processes. 
By  taking  them  into  account  one  obtains  a  greatly  improved  estimate  of  £ ,  at  the  ex¬ 
pense  of  computational  loads  depending  on  n.  On  the  other  hand,  if  the  variance  is 
also  ignored,  i.e.  the  multiplier  in  (2.  10)  is  deleted,  one  obtains  the  estimate  known 
in  numerical  analysis  as  Prony’s  method.  It  is  inconsistent  for  any  finite  A. 

From  (2.  7)  we  have  the  identity 


y. 


2lVT0 


=  d. 
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between  the  parameters.  We  can  therefore  estimate  ( 0 ^  0 by 


2|<S1  +  S() 


(2.  11) 


We  still  get  consistency  if  we  replace  Sh  by  Ah  in  and  by  2Bh  in  0*.  At  first  sight 
it  might  be  expected  that  this  alteration  would  lead  to  better  estimates,  because 
would  be  independent  of  the  remaining  random  variables  appearing  in  0*.  This  is  not 
the  case.  In  fact,  the  variances  are  uniformly  larger  than  those  of  the  estimates  as 

t 

written. 

Method  III  by  itself  cannot  distinguish  between  the  interval  (0,  n)  and  (it,  2tt)  of 
angular  frequency.  However,  by  combining  (2.  9)  with  the  first  cross -covariance  we 
can  consistently  estimate  co  (mod  2ir).  In  general,  we  have  for  any  fixed  integer  h 


n-h 

i  I\Vh 


t=l 


sin  2  <p6. jCOS  hw  if  -  0,  ±7r,  ±27t,  . .  . 

\  (2.  12) 

0  sinhw  otherwise  . 

s.  ^ 


Thus 


GO*  =  < 


cos  in  (0, 7r) 


cos  *  +  n 


if  C,  a  0 
1 


if  C,  <  0 
1 


converges  strongly  tow  /  kTt.  The  sampling  properties  of  this  estimate  become  poor 

-1  2 

when  go  is  near  a  multiple  of  7r,  because  the  variance  of  cos  £*  is  l/sin“oj  times 


_ _  2 

tThe  altered  version  of  0*  was  previously  suggested.  Cf.  Q  (£)  on  p.  8  of  22 L- 7529  and 

P.  .  of  (8.  3)  below. 
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that  of  £*,  and  the  latter  does  not  reach  0  at  these  points. 

If  cu  is  not  a  multiple  of  n,  we  can  combine  (2.  7)  and  the  second  line  of  (2.  12)  to 
give  another  pair  of  consistent  estimates: 


IV. 

"UV'’.* 

»*=  ■ 

J  S2  +  (2C1)2 

This  method  is  preferable  to  Method  III  when  A  is  sufficiently  small  (depending  on  £  ), 
but  the  relative  asymptotic  efficiency  approaches  0  as  A  increases.  If  a;  is  a 
multiple  of  it,  then 

/  2  ' 

0*  —  v  1  +  sin  2  (p 

as  n  — ■  06 .  Thus,  unless  the  unknown  initial  phase  in  (1.  2)  is  also  some  multiple  of  n, 
the  estimate  is  inconsistent. 

3.  AMPLITUDE  MOMENTS 

Here  we  derive  certain  formulae  which  will  be  needed  to  compute  the  large 

t 

sample  distributions  of  Methods  I  and  II.  The  reader  who  is  already  familiar  with 
the  properties  of  the  Rice  distribution  and  its  moments  can  skip  to  Sec.  4,  and  refer 
back  to  this  section  when  it  is  indicated. 

It  is  convenient  to  work  with  the  parameters  (0 A)  when  doing  the  analysis,  and 
to  simplify  the  writing  we  will  usually  drop  the  subscript  1.  Thus, 

0=2a2  A  =  p2/2<j2.  (3.1) 

Throughout  we  adopt  the  functional  definitions  used  in  [  1  ]  .  Any  symbol  [x.  y.  z.  ] 
refers  to  the  equation  so  numbered  in  this  handbook. 

t  Although  most  of  these  appear  in  existing  literature,  we  include  the  following  derivations 
in  order  to  have  a  permanent  record  of  verified  results. 
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Let  x  and  y  be  any  independent  normal  random  variables  with  common  vari- 
2  1 

ance  a  ,  and  respective  means  and  for  which 


2  2  2 

“it  +  ,i2t  =P 


for  all  integers  t.  The  joint  probability  element  of  (x^,  y  )  is 


-  -^j[(x-p1tr  +  (y-M2t)Z  ] 

^  e  dxdy. 


277(7 


2 


Changing  to  polar  coordinates,  the  probability  that  r  fall  in(r,  r  +  dr)  and  p^  in 
(p,  p  +  dp)  is 


-A  u,  cosp  +  u,,  sinp 

e  0  It  K  2t  K  , 
dF  =  — -re  e  drdp 

77  9 


2g. 


u.  = 


it 


it  9 


r  (i  =  1,  2). 


Writing  =  Rcosa  and  =  R  sin ot,  periodicity  gives 

,27 r 


,277  cos  p  +  u^sin p 


I  • 1 

J0 


dp 


-  f  pR 

-O' 


cos  c p 


d  cp 


[9.6.  16] 


p27 r 

■  v 


R  COS  ^dcp  =  277  I0(R). 


The  series  expansion  for  the  modified  Bessel  function  of  the  first  kind  of  general  order 
v  is 

V  (  /o\2n  +  l' 

[>•‘•“1  Vz)  •  Z  nirinVv"  T)  •  <3-2) 

n=0 


in 


Thus,  after  integrating  the  phase  out  of  dF  we  obtain 


f(r;e.*>  =^f—  re  *  i0(2vfr) 


(r  ^  0) 


(3. 3) 


I  2  2 

for  the  time -independent  probability  density  function  of  r^  =  v  +  y^  . 

If  we  substitute  I^'s  infinite  series  we  find  (3.3)  expressed  as  a  "mixture”  of 
central  chi  densities  with  Poisson  weights: 


wherein 


f(r;0,  \)=  £  cn(A)f2n+2<r;«) 


(3.4) 


n=0 


...  -4  x" 

C„W  =  e  JTT 


Vr;9>=' 


rt^>  e2 


V-l  0 
r  e 


The  latter  is  the  probability  density  function  of  n/  0/2  •  It  might  be  remarked 

in  passing  that  the  form  (3.  4)  generalizes  to  many  dimensions:  The  density  of 

I  2  2s 

r  =  w  x,  +  •  •  •  +  x  ,  when  the  x's  are  independent  and  x.  is  normal  with  mean  u. 

1  2  q  1  1  2 
and  variance  a  ,  is  given  by  (3.  4)  with  the  subscript  2n+2  replaced  by  2n+q  and  p~ 

2  2 
redefined  as  u  .  +  •••+  u  . 

th  q 

The  kLn  moment  of  y  is 

v 


—  y 

TT  -1  J0 


k 

2 


1  2 

2  ‘  v  '°  r(f) 

22  r(f)  2 
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Hence,  with  the  parameter  dependence  not  shown 


cv 


OC  — 

=  C  rkf(r)dr  =  ( 0)“  ^  c 

Jn  L,  n 


c  k 


r  x, 


2n+2 


n=0 


7  .  \  H(n+  7  +  I) 

=  fl2e'AX  - 2. - A" 

Z-/ 


(3.  5) 


n=n 


(n!) 


For  theoretical  purposes  it  is  best  to  express  O'  in  terms  of  Kummer  functions 

K 


[  13.  1.  2] 


M(a,  b,  z) 


V  r(n+a) 

L,  r(a) 


n=0 


F(b) 

T(  n  +  b) 


n 

z 

n  I 


Thus  k 

"k  =  °2  r^|  +  1)e"A  M(f  +  1>  *.*) 

k 

[13.1.27]  =  02  T(|  +  l)  M(  -  7,  1,  -A). 


(3.  6) 


The  latter  form  agrees  with  Rice's  formula  (p.  339  of  Wax  [  7  ]),  who  uses  the  confluent 
hypergeometric  notation  Fj  in  place  of  M. 

The  formula  for  the  even  moments  reduces  to  a  Laguerre  polynomial.  In  fact, 


[13.6.9] 


o?2k  =  0kk!  M(-k,  1,  -A)  =  0kk!  Lk(-A) 
k 

Lk<-«  ■  I 
1=0 


r k  \  a^ 

\  j  /  i  ! 


(3.7) 


This  can  be  derived  directly  from  (3.  5)  by  first  noting  that 
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and  then  using  Rodrigue's  formula. 

We  examine  the  behavior  of  the  amplitude  moments  for  large  A.  From  (3.  7)  we 
immediately  see  that 

“2k  s  <M)k 


as  A  — ■  oc.  For  the  odd  moments  we  use  the  expansion 


k-.1 

[13.5.1]  M(-|,  1,-A) 


N-l 


iril^nFir  +  o<*  "> 

n=0 


,2  A 


-n 


y 

j 


wherein 


r(n+|- 

r(i  -k) 


Since 


[  6.  1.  17] 


1 

r(k+|)  r(*  -k) 


i 

—  sin  7t(k+ 
7 r 


we  find  for  (3.6),  after  replacing  k  by  2k-l, 


2k-l 


rz(k+^) 


°2k-l  “  (6X) 


1  + 


N-l  2 

r  (n+i-k)  _ 

-  A  “  +  0(A 

7tn  ( 

n=l 


n  . 


)  (3.8) 


th 

as  A  —  °°.  Consequently,  the  k  moment  of  r/p  tends  to  1  as  A  increases,  which  is 
obvious  for  9  tending  to  0  with  p  fixed.  In  a  sense  made  precise  by  Rice  on  p.  240  in 
Wax  [  7  ]  ,  the  density  of  r  looks  normal  about  p  with  variance  \  9  (not  9)  when  p  is 
large  and  9  is  fixed. 


We  will  need  <v  and  tv  ,  made  independent  of  0  after  diversion  by  o  ^  raised 
to  the  appropriate  power.  For  computational  purposes,  it  is  best  to  write  these  as 
linear  combinations  of  I()  and  l  ^  The  connection  is  made  via 

A 

M(  2,  1,  A)  =  e2  IQ(|) 

A 

[13.6.3]  M(|,3,A)  =  e2 

and  certain  recurrence  relations  among  3  M's  at  a  fixed  argument.  From  the  first 

formula  in  (3.6)  we  see  that  ot  ^  requires  M(-y,  1,  A)  and  requires  M(tj,  1,A).  If 

we  multiply  [  13.  4.  2]  with  a  =3/2,  b  =  2  and  z  =  Aby  A,  and  add  to  it  -2(1 +A)  times 

3 

[  13.4.4]  with  a  =3/2,  b=l  and  z=A,  we  obtain  a  formula  for  M(— ,  1,A)  in  terms  of 
the  above  two  M's.  The  result  is 

A 

Q?1  =  e  2  [(1  +  A)I0(|)  +  AI1(|)].  (3.9) 

5  13 

We  get  M(  — ,  1,A)  in  terms  of  M(— ,  1,  A)  and  M(  — ,  1,  A)  directly  from  j  13.4.  1]  with 
a  =  3/2,  b=l  and  z=A.  Thus, 

3  _A 

«3  =62  ~L  e  2  [(1  +  2A+|a2)I0(|)+|a(2  +  A)I1(|)];  (3.10) 

or,  ot  can  be  written  in  terms  of  I(^  and  a 
4.  ANALYSIS  OF  METHODS  I  AND  II 
From  (3.  7)  the  mean  power  is 

a  =  0(1  + A) 

and  thus  by  (3.  9) 
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O' 


F(X)  = 


O' 


\ln 


X 

2 


(1  +  X) 


[(1  +  X)I0(|)  +  XI1(|)] 


Using  (3.  2)  and  (3.  8)  we  see  that 


(4.1) 


F(0)  =  n/tt/2  =  .  8862  2693 


F(oe)  =  1. 


Since 

[9.6.28] 


d^Vz>  =  ‘i<z) 


il1(2)  =  I0(Z)-|l1(z) 


we  find  for  the  derivative  of  F 


^ir 


X 
'  2 


F'(X)  =V 


4  (i+x)3/2  12 


Mt)>o 


(4.2) 


Thus,  F  is  a  strictly  increasing  function  of  its  argument,  and  the  estimate  (2.  3)  of  X 
is  indeed  well-defined. 

To  compute  the  large  sample  distribution  of  (2.4)  we  need  that  of  (a^,  a ^),  of 
which  (6*.  6*  )  are  functions.  In  general,  it  is  true  that  (Cramer  [  2  ] ,  p.  364) 


''"‘VV-  -^VV . '/n<‘‘k  '"k  * 

11  2  2  p  p 

tend  to  joint  normality  as  n  — ■  00  with  zero  mean  vector  and  covariance  matrix  with 
..th 

ij  entry 


k.+k. 
i  J 


a;  a. 


k. 

J 


“k,1  <V 

1  ] 


",  ) 


(i,  j  -  1,  3, . . . ,  p) 


(4.3) 
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In  particular,  for  vn(a  -o'  )  and  \n(a  -o'  )  we  have  the  covariance  matrix 

11  2  2 


A  = 


2 

O'  -  O' 

2  1 


0?  -  O'  O' 

3  1  2 


O'  -  01  01 

3  1  2 


O'  -  ct 
4  2 


We  express  the  entries  in  terms  of  F(A)  and  (see  (3.  10)) 


a 


G(A)  = 


o? 


3/2 


3\/' 


7t  e 


A 

2 


4  (i+X)572^+2-f^2)1°(l)+>^>Ii<l)1 


We  then  have,  using  (3.  7), 


2 

o;  -  ct 
2  1 


=  0(1 + A) (1  -  F  ) 


a  -a  a  = 
3  1  2 


2 

a;  -  o; 

4  2 


03/2(i+a)3/2(g-f) 


=  0  (1  +  2A) 


where  the  argument  value  is  understood. 

We  apply  the  delta  method  (see  Appendix  A),  and  introduce  the  functions 


1  +  F  ^Xj/x*) 


g  =  g<x,.x2)=( 


if  x^/xj^  —  ^ln/2 


otherwise 


Si  =  g 


&2  X  2 " ^ 


(4.4) 
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In  this  notation  the  estimates  (2.  4)  are 


S1  =gl(al'a2> 


°2  =  s2<ara2) ' 


Under  the  assumption  X  f  0,  i.e.  » 1/«“  >  n/tt/2,  the  matrix  of  derivatives  evaluated 
at  the  true  parameter  values  is 


D  = 


9gl 

9xl 

9  g2 
9X1 

II 

hi 

-hl 

ffl 

3x 

*2 

9x2 

h2 

l-h2 

_  z 

z  __ 

X  =  a? 

L  — 1 

wherein 


h  =  -2£. 
1  8x 

1 


x  =  a 


-e 


(l  +  \)3/2F' 


h  = 

2  3x„ 


X  =  O' 


+ 


1+x  2(1  +X)2  F' 


The  covariance  matrix  of  the  limiting  normal  distribution  of  \ln(0*  -  dj  and  N ln(0 *  -  Oj 


is 


2  =  D'  A  D  =  0  Q 


(4.5) 


where  the  symmetric  matrix  Q  -  [Q  ]depends  only  on  X.  We  express  its  elements  in 
terms  of  F, G  and 


H  =  (1  +  X)F' 


(4.6) 
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After  considerable  algebraic  manipulation  we  arrive  at 


The  determinant  of  Q  is  easy  to  compute  because  that  of  D  is  simply  h^.  We 


have 


2  2/2 

(1  +  A)(1-F  )  (1  +  A)  '  (G-F) 


(1  +  A)'V2(G-F) 


1  +  2A 


[(1  +  A)(1  +  2A)(1-F2)  -  (1  +  A)3  (G-F)2]  . 


(1  +A)H 


Consequently, 


(4.  8) 


The  forms  (4.  7)  and  (4.  8)  are  suitable  for  analytical  examination  for  small  and  large 
argument  values,  which  is  done  in  the  next  section.  The  problem  of  numerical  evalua¬ 
tion  is  discussed  in  Sec.  10. 
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The  consistency  of  Method  II’s  0^  -  estimate  in  (2.5)  follows  immediately  from 
the  formula  for  the  fourth  moment  of  r  given  in  (3.  7).  Indeed, 

«4  =  202(1  +  2A  +  ^A2)  =  202  +  40l  02  +  02  =  2«2  -  02 


because  0 .  +  0„  is  o>  . 

12  2 

To  compute  the  large  sample  distribution  we  employ  precisely  the  same  tech¬ 
nique  used  for  Method  I.  Since  the  estimates  now  depend  on  (a?,a4),  rather  than 
(a^.a^),  we  have  a  much  simpler  A  matrix.  From  (4.3)  and  (3.7)  we  have 


a . 


a 


a  -  a  a 
6  2 


2 

a  -  a  a  a  -  a 
6  2  4  8 


with 

a4  -  a2  =  02(1  +  2A) 

3  2 

a  -  a  a  =  40'  ( 1  +  3A  +  A  ) 

6  z  4 

a4  =  404(5+2OX  +  13x2  +  2x3)  • 


For  the  g  function  we  take 


g  =g(x2,x4)  =  ^ 


(2x2  "  X4>2  if2x2~X4 


otherwise  , 


so  the  estimates  (2.  5)  are  0*  =  &2  -  g(a2>  a4)  and  0*  =  g(a2>  a4).  Again  assuming 
A  ^  0,  the  matrix  of  partials  evaluated  at  the  true  parameter  point  is 
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D  = 


1_h2  h2 


-h.  h, 
4  4 


wherein 


h  = 

2  ax2 


X  =  OL 


- 2(1+^ 


h  = 

4  9x 


x  =  O' 


20  X 


The  form  (4.  5)  again  holds.  The  elements  of  Q,  after  many  frustrating  polynomial 
manipulations,  are  found  to  be  simply 


(4.  10) 


The  determinant  is 


(4.  11) 
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The  two  variances  and  the  squared  correlation  are  graphed  in  Figs.  2,3  and  4.  The 
corresponding  curves  for  Method  I  have  the  same  qualitative  properties. 


5.  COMPARISON  OF  METHODS  I  AND  II 

In  this  section  we  compare  the  covariance  matrices  (4.  7)  and  (4.  10)  for  small 

and  large  values  of  the  parameter  A.  The  figure  of  comparison  is  the  ratio  of  their 

determinants  already  introduced  in  (2.6).  The  determinant  of  a  covariance  matrix 

is  called  the  generalized  variance  of  the  distribution  (centered  on  the  origin).  The 

reason  for  this  terminology  is  that  the  area  of  the  region  interior  to  the  ellipse 
-1  2 

x'  Q  x  =  c  ,  a  contour  of  constant  probability  in  the  bivariate  normal  case,  is 


7r c “  I  Q  1 2  (cf.  Cramer  [  2  ]  ,  p.  300). 

Let  us  examine  for  small  values  of  A  first.  From  (3.  2) 


I^(z)  =?  (z/ 2)1  /T(v  +  1)  as  z  —  0. 


Referring  back  to  (4.  1),  (4.  2),  (4.  4)  and  (4.  6)  we  have 


F(0)  =  ^7r/2  G(0)  =  3^tt/4 


(5.1) 


H(A)  ^  F'(A)  ^  n/ttA/16. 


From  (4.  7)  we  see  that 


22,1 


have  equal  limiting  values  as  A  —  0;  viz. 
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We  get  the  same  thing  for  H~  times  (4.8),  so  that 

,  ,  ih(1h-57T) 

-1  “  X2 

7TA 

Combining  this  with  (4.  11)  gives  the  value 

d(0)  =  1.49 


(5.  2) 


(5.  3) 


of  (2.  6)  at  the  origin. 


To  approximate  (4.  7)  and  (4.  8)  when  A  is  large  requires  a  bit  more  effort.  We 
use  the  formula  (3.  8)  with  k  =  1  and  2,  and  keep  terms  up  through  1/A'  .  Setting 
e  =  1/A  we  have 


F(A) 


=  (1  +  e)'if  1  +  F  (30/2)  [  r2(3,)e  +  *  T2(3/2)e2  +  -  T2( 5/2)e3  +  0(e4)] 1 


l 


7T 


“i-»€+He2'l§f3  +  0<£’) 


(5.  4) 


and 


G(A) 


=  (1  +  e)  2  /i+r~(5/2)[r2(  ->)c  +  ir2(-*)e2  +^r2(3/2)e3  +  0<e4)]) 


.  3  39  2  209  3  .  4 

=  1+^-He  +  +0(€  } 


as  e  — ■  0.  It  is  not  true  that  the  asymptotic  expansion  of  F'(A)  is  given  by  the  result 
of  differentiating  (5.  4).  We  must  return  to  (4.  2)  and  use  the  large  argument  expansion 
[  9.  7.  1]  for  1^.  In  terms  of  e  we  obtain  for  (4.  6) 


4H  =  €  [  1  -  ^ e  +  ~  e2  +  0(e3)]  . 


Thus, 


22 


2(1-F“)  =  e  [  1  e  + 1  e2  +  0(e'3)] 

G-F  =e[l-|e+|2e2  +  0(e3)] 

2(H  +  |F)  =  1  +  U  --^e2  +  0(e3)  . 

After  some  series  manipulations,  we  find  the  quantities  in  (4.  7)  are  in  the  limit  the 
same  as  the  corresponding  ones  in  (4.  10): 

Q.. 

lim  =  1  d(°°)  =  1  .  (5. 5) 

\  gij,II 

2 

If  we  retain  terms  up  through  e  in  the  square -bracketed  quantity  in  (4.  8),  we  find 
in  terms  of  A 


|  Q  |  =*  0(1/X)  +  8  +  4\  (5.6) 

as  \  —  00 .  This  is  to  be  compared  with  the  linear  approximation  9  +  4X  to  |  |  . 

The  determinants  (4.  8)  and  (4.  11),  and  their  ratio,  are  plotted  in  Figs.  5  and  6. 

6.  ASYMPTOTIC  EFFICIENCY  USING  AMPLITUDE  DATA 

Let  (0*,0*)  be  estimates  of  (flj,#^)  basec*  on  independent  observations 

r  .,  r_, . . . ,  r  having  a  common  probability  density  function  f,  and  suppose  \ln(6*-6  ) 

1  z  n  11 

and  N/n(O*-02)  tend  to  0  mean  normality  as  n-*  with  covariance  matrix  2  .  The 
joint  asymptotic  efficiency  of  (0*,0*)  is  defined  as 


(Cramer -Rao)  where 


e  =  ■ 


L  I  Is 


<  1 


1=- 


t 


9  2  logf 


de.dd. 

i  j 


(6.  1) 


(6.2) 
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is  the  so-called  information  matrix.  The  expectation  is  taken  under  the  assumption 

t 

that  is  the  true  parameter  point,  and  both  J_  and  2  are  presumed  nonsingular. 

The  regularity  conditions  which  ensure  the  existence  of  a  consistent  solution 

A  A 

(0^,  0,;)  of  the  Likelihood  equations 


9  lQg  L  =  n  al°sL  =  0 

9 6  9 0 

1  2 


n 

L*TI«VW 

t=l 


such  that  \/n(0^ -0^ )  andN/n(09  -0,f)  are  asymptotically 

N(o,  2'1) 

are  satisfied  by  the  density  (3.  3).  (We  indifferently  write  f(  •  ;  0^,  0.})  and  f(  •  ;  0,  A) 
for  the  same  function. ) 

It  is  important  not  to  infer  any  optimum  fixed  sample  size  properties  of  the 
Maximum  Likelihood  (abbrev.  ML)  estimate.  Indeed,  since 

Io(2/lr)i'Hl(/l)  H2«r)’ 

there  exists  no  jointly  sufficient,  and  hence  no  jointly  efficient,  estimate  of  (0 ^ ,  0  .,) 
(except  when  A  =  0).  The  practical  implication  of  this  is  that  we  are  not  guaranteed 
a  unique  pair  of  solutions  of  the  Likelihood  equations  for  any  fixed  n.  However,  as  we 
shall  see,  there  is  a  unique  pair  of  nontrivial  solutions. 

We  continue  to  write  the  various  formulae  in  terms  of 

9  -  q  x  =  vci 

t  For  the  general  material  used  in  this  section  see  Sec.  32.  6  of  Cramer  [  2  ]  and 
Secs.  18.  22  -  18.  28  of  Kendall  and  Stuart  [  4  ] . 
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and  introduce  the  abbreviation 


R  =  r 


',(*71' 


IJ2.  Tr 


(6.  3) 


From  the  logarithm  of  (3.  3)  we  then  have 


6  al°Sf  =  -2  ;R+f+H 

30,  \/  0  0 

1 


iissl  =  — - —  r  -  i 


30, 


2  \/0A 


(6.4) 


The  Likelihood  equations  are 


o  =  = -2  jl  )  r 

30.  V  0  Z_y  t 


na 

+  _£  +  n(A-l) 


t=l 


n 


0  =  <,25i£k  =  V" 

Ml  -Jex  L‘,  - 

t  =  l 


th 


wherein  R  means  R  with  r  replaced  by  the  t  observed  value  r^.  If  we  multiply 
the  second  of  these  by  2A  and  add  the  result  to  the  first  we  obtain  0(1  + A)  =  a.)t  or 


VS2  =  a2 


(6.5) 


in  terms  of  the  quantities  of  interest.  The  other  Likelihood  equation  is 


n 


T<X)=- 


2r  x 
t 


a2'X 


t=l  L 


2r  x 
t 


=  x. 


(6.6) 


a2  "  X 


V 
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Since  I  vanishes  at  the  origin,  x=0  is  always  a  solution  (which  is  incorrect  unless 
1  2 

A  =  0).  There  is  one  other  root  x  =  6^  which,  together  with  (6.5),  gives  a  consistent 
asymptotically  efficient  pair  of  estimates  of  (0^,0^).  We  can  approximate  the  desired 
solution  by  successive  substitutions 


x 


i+1 


T(x.), 

l 


initialized  by  (say) 


t 

with  A*  given  in  (2.3). 

c  k 

We  first  show  that  (  R  exists  for  every  fixed  k,  and  in  the  process  of  so  doing 
obtain  the  explicit  formula  for  (  R.  After  changing  the  dummy  of  integration  to 


z  =  2 


we  have  integrals  depending  on  A, 


k+2 

=  (jA)_  V 
Yk  k 


a  = 
k 
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2 


I()(z)dz 


(6.7) 


where  has  already  been  calculated  in  Sec.  3.  The  ratio  I^(z)/Ig(z)  increases 

t  This  procedure  might  run  into  difficulty  when  the  underlying  parameter  A  is  small. 
A  different  method  for  approximating  the  nonzero  root  of  (6.  6),  which  circumvents 
this  problem,  is  used  by.T.  R.  Benedict  and  T.  T.  Soong  of  Cornell  Aeronautical 
Laboratory. 
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steadily  from  0  to  1  with  z.  Thus,  since  the  integrands  are  positive  functions, 

2 


7.  s  zk+1e  4A  I  (z)dz 


(6.8) 


for  all  k  l.  Since  IJ(z)dz  =  dl^(z),  and 


[9.7.1] 


Vz)  *  — 

N  2nz 


as  z 


for  all  orders  v,  we  find  in  the  notation  (6.  7) 

2 


oc 

£ 


k+1 
z  e 


z 

4A 


y, 


1 1(z)dz 


k+1 

2A 


-  (k+l)y  ,  <  00 
k-1 


(6.9) 


Thus,  all  moments  of  R  exist.  For  k=l  there  is  equality  in  (6.8),  and  consequently 


<5r  = 


20 


,  .3/2  A 

(4A)  e 


2A  T0 


=  ^0  A 


i.  e.  p 


(6.  10) 


because 


%  =  1  01 2  =  ^4  +  ^' 


In  the  same  way  the  change  of  variable  and  integration  by  parts  gives 

r2 

2.  ~  -X  »“ 


c  2  2e"  f"4  "»,/',/» 

frR-— Ve  ’A2'/?1 


dr 


20 


3/2 


.5/2 


y. 


(4A) 

=  0  ^0  A  (2  +  A) 


2T  -4y2 


(6.  H) 
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because 


a  =  262(1  +  2X  +  ‘  X2). 


From  (6.  4)  and  (6.  10)  we  get 


C  9  1<>K  f  =  f  9  loS  f  =  0 

c  do  c  do 


d  0.. 


From  this  it  follows  that  (6.  2)  is  alternatively  given  by 


1  = 


f  91°gf  ai°g 

C  Xa  P)fl 


do.  do. 

1  J  J 


(6.  12) 


which  is  easier  to  work  with  for  the  particular  density  under  consideration  (although 
in  most  problems  (6.  2)  is).  Using  (6.  10)  and  (6.  11),  together  with  the  by  now 
familiar  formulae  for  a  and  a  ,  we  obtain  from  (6.  4)  and  (6.  12) 


<9ZJ  =  —  f  R2  +  1  -  2X  -4X“  =  1  -  2X  +  4X^4> 

J11  0  l 


e2j  =  -£/V  +  i  +  2x 
J 12  0 


=  1  -  2X4> 


fl2  =  J_^R2-1 
6  J22  ex1  R 


=  <f>  . 


(6.  13) 


We  have  put  (see  (6.  7)) 


<f>  +  1  = 


ex 


_£*  r~  3  'lx  'i(z)  , 

8  X3  0  Z  '0<Z>CZ 


(6.  14) 


c2  c  2 

which  always  exceeds  1  because  0X  =  q  R  —  £ R  . 
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Since  4>  is  independent  of  0  we  can  write,  in  keeping  with  our  previous  notation. 


2?  =  J__1  =  02Q° 


under  the  assumption  that 

- L_  =  e4  |j|  =  (1  +  2X)  $-1  t  0.  (6,15) 

I  Q°  I 

We  then  have  from  (6.  13) 


ML  (amplitude). 

<i> 

2X*-1  n 

D°  ^ 

^  (1  +  2X)4>-1 

2X  4>-l 

1-2X  +  4X“$] 

for  the  covariance  matrix  of  the  limiting  normal  distribution  of  ^n(  6^- 0^)/0^  anc^ 

v'n(0  -0  1/0-  The  function  $=  4>(X)  defined  in  (6.  14)  decreases  steadily  from  1  to 
2  2'  1 

0  as  X  goes  from  0  to  Figure  7  graphs  the  results  of  doing  the  integral  numeri¬ 
cally,  and  Fig.  8  the  ensuing  values  of  the  determinant  of  Q°.  See  Appendix  B  and 
the  acknowledgment  contained  therein. 

We  can  check  the  calculations  leading  to  Q°  by  using  (6.  5).  If  we  square 


x/n(0  -0^  +  n/ n(02-02)  =  ^n(a7  -  ) 


and  take  expectations,  we  have  (see  (4.3)) 


2Q12  +  Q22 


1  +  2X 


in  the  limit  as  n  —  oc.  Upon  substituting  the  formulae  for  the  Q°  ’s  we  obtain  the 
desired  redundancy,  =  4>. 
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We  can  deduce  expressions  for  the  elements  in  (6.  16)  when  A  is  large  by  the 
following  trick.  It  was  pointed  out  following  (3.  8)  that 

r  is  approximately  N(p,  3,0) 

for  r  in  the  neighborhood  of  large  p  with  0  fixed.  Let  us  suppose  we  are  actually 

sampling  from  this  distribution.  The  ML  estimate  of  3,0,  when  the  mean  is  unknown, 

is  the  sample  variance  =  a^  -  a”  .  As  n  — •  00 ,  \l n(m^  -  p.;)  tends  to  be  normally 

distributed  about  0  with  variance  p  -p“  which,  by  the  underlying  normality,  is 
2  .  4 

equal  to  2p.).  Since  p^  is  £0,  this  is  tantamount  to  saying 

\l~n(0  -0)  ~  N(0,  20  “) 

where  0,  as  always,  is  0^  .  Referring  back  to  (6.  16),  we  should  therefore  have 


j> 

(1  +  2A)  tf-l 


2. 


Solving  for  $  gives 


Thus 


<f>  3 


2 

1  +  4A 


_1_ 

2A 


-  +  o(a'3). 

8A 


-1 

1  +  2A 


as  A  — ■  00 , 


with  linearly  increasing  determinant 

|  Q°  |  &  1  +  4A  . 


(6.  17) 


(6.  18) 


(6.  19) 
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It  is  shown  in  Appendix  B  that  (6.  17)  is  indeed  the  case,  and  consequently  both 

Methods  I  and  II  will  approach  full  asymptotic  efficiency  as  X  increases.  If  we 

2 

choose  2  as  2^  =  0  in  (6.  1)  we  obtain  e,  as  expected,  as  a  function  only  of  X: 


enw  = 


kj°(x)  I 
lQn(x)| 


(6.  20) 


We  write  the  asymptotic  efficiency  of  Method  I  as 


en(X) 

ei(A)  "  “d(xT 


(6.21) 


where  the  relative  measure  d(X)  has  already  been  shown  in  (5.  5)  to  equal  1  at  X  -  00 
Thus,  upon  comparing  (6.  19)  and  (4.  11), 


ej(«°)  =  ejj(°°)  =  *• 


We  note  in  passing  that  the  estimate  of  6  given  by  Method  I  also  becomes  the 
sample  variance  when  we  replace  F  by  the  large  argument  approximation 


F(X)  = 


given  by  (5.  4).  Inverting  this,  we  have  from  (2.  3)  and  (2.  4) 


0* 

1 


1  + 


(ai/a2> 

l-(a~/a2) 


a2  ”  ai 


For  small  values  of  the  parameter  we  have  (see  Appendix  B) 


(6.22) 
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=  1  -  2A  +  5A2  +  0(A3) 


(6.  23) 


which  combines  with  (6.  15)  to  give 

lim  A2  |  Q°  |  =1.  (6.  24) 

A  -0 

Thus,  from  (4.  11)  and  (5.3) 


e(0)=.68  e  (0)  =  1.  (6.25) 

The  functions  (6.  20)  and  (6.  21)  are  graphed  in  Fig.  9. 

7.  DISTRIBUTION  OF  SAMPLE  COVARIANCES  AND  CROSS  COVARIANCES 

We  now  turn  to  somewhat  different  theoretical  considerations,  viz.  the  joint 
large  sample  distribution  of  the  random  variables  appearing  in  (2.  7)  and  (2.  12).  In 
anticipation  of  future  needs,  we  consider  a  pair  of  more  general  time  series 


x 

t 


F  +  u . 
t  It 


y  =  G  +  u_ 
yt  t  2t 


(7.1) 


where,  for  notational  simplicity,  we  allow  t  to  range  over  all  positive  and  negative 
integers.  We  take  {u  }  and  {u  }  to  be  independent  and  identically  distributed  as 

-L  L 

some  zero  mean  linear  Gaussian  process,  say  {u^},  with  a  covariance  sequence 


00 


h=-oc 


(7.2) 


The  following  assertions,  for  the  class  of  "uniformly  asymptotically  stationary" 
regression  sequences,  are  proved  in  Appendix  C. 
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Theorem.  Let  {F^}  and  { G^}  in  (7.  1) 


be  any  two  bounded  real  number  sequences 


(t  =  0,  ±1,  ±2,  . . . )  such  that 


n 

-N  F  F  =  x[F)  +  0(1/ N/n)  (7.3a) 

n  Li  t  t+h  h 

t=l 

n 

-  \  G  G  =  xjpV  o(  1/N/n )  (7.3b) 

nL  t  t+h  h 

t=l 

n 

-  \  F  G  ,  =  /i.+ o(l/N/n)  (7.3c) 

n  ^  t  t+h  h 

t=l 


where  the  order  terms  exist  independently  of  h,  and  the  X* s  are  necessarily  even. 
Define 


n 


h  n/_j  t  t+h 
t=l 


r(n)  =  1  >  x  v 
h  n/z  t  ^t+h  ‘ 
t=l 


(7.4) 


Then,  for  every  fixed  integer  h, 

,(n) 


A'  '  -* 


A(F)  +  a 
h  h 


<n)-  % 


(7.  5) 


as  n  —  00  both  in  mean  square  and  with  probability  one.  Set 


lZn)  =  \/n(A<n)-  A(F)-  a  )  v/n)=  -s/n (C<n)-p  ) 

h  h  h  h  h  h  h 


Then 


Urn  -  Z  ffk(2<A<k+b-a  +  4tL)+  <Vb-a  + Vb+a>] 


(7.  6a) 


k=-°o 
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lim  l  W(n)W(n)  =  N  a  [a£>  +  >K  +(T.+,  1 

n  (  a  b  i_j  k  k+b-a  k+b-a  k+b-a 


k=-°° 


lim/  l/1^  w|n^  =  ^  a [g  +Mlr,w,  1 

„  a  b  /_,  kl  k+b-a  k+b+a  J 


(7.  6b) 


(7.  6c) 


k=-° 


.(F) 


for  any  integers  a  and  b,  provided  2  a  A  etc.  exist.  Furthermore,  any  collection 

K.  K. 

of  random  variables 

rj(n)  i/n^ 

Uh/-",Uh  ’  Wh  ’  ’  ’  *  ’  n 
1  p  p+1  p+q 

have  a  large  sample  (p  +  q)  -variate  normal  distribution  with  zero  mean  vector  and 
covariance  matrix  whose  entries  can  be  obtained  from  the  formulae  (7.6). 

We  further  set 


v*n)  =  "MBJ,n)-  \G)  - )  Bin)=  { )  ytyt+h  <7-  7> 

h  h  h  n  h  n/_,  t  t+h 

No  separate  proof  is  required  to  establish  that 


t=l 


lim  /  V(nVn)  =  (7.  6a)  with  (F)  replaced  by  (G)  (7.  8) 

^  a  b 


and 


lim  f  U(l,)V<n)  =  0 
a  b 


(7.9) 


,(F). 


by  x  and  y  independence.  Since  (7.6c)  does  not  depend  on  the  A  s  which  appear  in 
the  centering  of  the  U's,  it  follows  that 

lim  £v(n)  W(n)  =  lim  U(n)  W^n)  (7.10) 

n  a  b  a  b 

for  all  a,b(  as  can  be  verified  by  direct  evaluation).  V's  can  be  added  to  the  statement 
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concerning  joint  asymptotic  normality  of  U's  and  W’s,  and  the  preceding  formulae  used 
to  compute  covariances  of  the  distribution. 

Specialization.  Our  model  (1.1)  is  subsumed  under  (7.  1)  with 


F  =  p  cos(wt  +  tp)  G t  =  P  sin(cat  +  cp) 


ch  =  *  6h,  0- 


(7.11) 


It  is  shown  in  Appendix  D  that  the  hypotheses  in  (7.  3)  hold,  for  any  h  <  0,  with 


(F)  2  2  (G)  .  2  2 

A,  =  cos  <p  •  p  coshoo  A  =  sm  cp •  p  cos  ha) 

h  h 

if  a)  =  0,  ±ir,  ±2^  . . . 


3F>=*<G>=  ip2cosha, 
h  h  2 

if  u>  f  0,±7T,±27T, . . . 


(7.  12) 


and 


Mh  = 


=  < 


|  sin  2 cp’  p  cosha) 


1  2  .  . 

2  p  sin ha> 


if  o>  =0,  ±7r,  ±2tt, 


if  w  f  0,±7r,±27r, 


(7.  13) 


From  our  Theorem,  the  limits  asserted  in  (2.  7)  and  (2.  12)  are  thus  seen  to  be  correct 
(the  change  of  the  upper  limit  from  n  to  n-h  in  the  definitions  of  A,  B  and  C  having  no 
effect  in  the  limit). 


We  now  put 


R<n)  =  Wn)-  v  )  =  U<n)  +  V<n) 
h  h  n  n  n 


T[n)  =  N/n(2cJn)-  2a  )  =  2w[n) 
h  h  n  n 


(7.  14) 
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in  the  notation  of  the  Theorem.  Any  collection  of  random  variables 


(n)  (n)  (n)  (n) 

h .  h  h  ,  .  h 

1  p  p+1  p+q 

tends  to  joint  normality.  Using  (7.6)  through  (7. 13),  we  find 


lim /R(nVn)  =  02U(<5  .+5  )  +  2\cosawCosbw] 

n  ^  a  b  a,  b  a, -b 


(all  w) 


lim  /T(n)  T^n)=  62  [  6  +  2Xcos(b-a)w] 

„  u  a  b  L  a,b 


(7.  15) 


and 


lim  (c  R(nMn*  =  <| 
^  a  b 


20  \  cos  aw  sin  bw 


2 


2 

We  see  that  the  limiting  value  of  (C  is 


W  f  0,  ±7T,  ±27T, 


26 sin  2(/>cos  aw  cos  bw  w  -  0,  ±7r,  ±27r, 


(7.  16) 


2  2 
0  (1  +  2A.)  =  , 

i.e.  the  variance  of  <s/n(a  -a  ).  This  provides  a  check,  since  S0  is  &2  and  yQ  is 


W“2' 


8.  ANALYSIS  AND  COMPARISON  OF  METHODS  III  AND  IV 

We  again  apply  the  "delta  method"  of  Appendix  A  to  obtain  the  covariance  matrix 
2  of  the  asymptotic  normal  distribution  of  N/n(0*  -0^)  and  'v/n (0*  -6^)  for  the  estimates 
(2.  11).  To  economize  on  space  we  drop  the  notational  distinction  between  variables 
and  true  parameter  values  made  in  Sec.  4.  For  Method  III  we  set 
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g  =  g(yQ>  y^)  =  r2  -  2^yi  +  r0 


4  = 


y2  +  Tv^y 


4y, 


where  (see  (2.  7)  -  (2.  9)  in  6,  X  notation 


70  =  e(l  +  X)  yi  =  0A.4  y2=ex(2|  -1)  . 


Then 


and 


e1  =  g  =  g1 


e;=*i<so-si’s2> 


°2  y0 " g  g2 


^-*2<S0'S1-S2> 


There  are  here  two  functions  of  three  variables.  The  matrix  of  partials  is  found  to 
depend  only  on  4  : 


D  = 


®*1 

% 

ag2n 

% 

1  +  2| 2 

r 

o 

agi 

dg2 

1 

-4| 

44 

dy\ 

l  +  2£2 

3gi 
_  ay2 

dg2 

dy2_ 

1 

-1 

(7.  15)  gives  the  covariance  matrix  of  the  joint  asymptotic  normal  distribution  of 


^n(SQ-y0)  ,  ^(S1-y1)  ,  \/n(S2-y2) 
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It  is 


1  +  2A 

2A£ 

2A(2£  2  - 1) 

2A| 

\  +  2A£ 2 

2\^(2r  -i) 

2 

2 

,  ,  2  ,.2 

2A(2£  -1) 

2A|(2|  -1) 

i  +  2A(2|  -1) 

After  some  algebra  we  find 

2  =  D'AD  =  0“P 


(8.  1) 


where  P  =  \  P  1  depends  on  f,  as  well  as  X,  but  not  on  0.  The  dependence  on  the 

~  ij 

former  is  via  the  value  of 


1  +  16g  2 

2  2 

(i+2  rr 


(8.2) 


and  the  entries  are  simply 


III. 

P11  = 

1  +  K 

P12  = 

-K 

P22  = 

2A+  K 

(8.3) 


The  generalized  variance  of  the  distribution  of  N/n(0*  -0^)/  0^  and  (0*  -  0^)/  0^ 
thus 

III. 

|  P  |  =  K  +  2(1  +  K)A 


is 


(8.4) 
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The  large  sample  distribution  of  the  cosine  estimate  (2.  9)  is  computed  by  the 
same  technique.  We  find 


N/n<r-  O  ~  N(  0,  — -  • 

'  2  A 

t  2 

without  restriction  on  ,  =  cosca.  X  times  this  variance  function  is  plotted  in  Fig.  10 
2 

against  ,  .  If  we  assume  ca  is  not  a  multiple  of  7r  and  invert  the  cosine  estimate  to  get 

2 

an  estimate  of  frequency,  the  variance  in  (8.  5)  is  to  be  divided  by  sin  ca. 

For  the  estimates  (2.  13)  we  obtain  a  covariance  matrix  which,  when  written  in  the 
form  (8.  3),  has 


i-3r  +  ^4  a 
<1  +  2, ¥  ) 


(8.5) 


IV. 

K  replaced  by  L=  1 +  2,  ~(1-,^)A 

¥  t  1) 


(8.  6) 


This  depends  on  X  as  well  as  ,  . 

2 

Figure  11  is  a  plot  of  (8.  2)  against  ,  .  It’s  minimum  value  at  0  is  \  and 
maximum  at  3/8  is  8/7.  Thus,  from  (8.4) 


min  I  -m  I  =  ^  +  3X 


max 

£ 


(8.7) 


It  is  shown  in  the  next  section  that  the  covariance  matrix  of  the  limiting  normal  distri¬ 
bution  of  ^n(  9^  -  9^)/ 9 1  and  (  ^  '  is  given  by  (8-  with  K  =  Thus, 

|  P° |  =  2A  .  (8.  8) 


t  The  variance  in  (8.  5)  is  5  the  one  given  by  Eq.  (7.  2)  in  [  3] .  This  is  as  should  be, 
since  there  we  had  available  only  a  single  time  series,  i.  e.  ,  *  was  computed  using 
A's  rather  than  S's. 
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Figure  12  is  a  graph  of  these  three  straight  lines. 

The  asymptotic  efficiency  of  Method  IV  relative  to  Method  III  is  less  than  or 
greater  than  unity  depending  upon  the  sign  of  the  difference 

|PIV|-  |  Pjn  |  =  (1  +  2A.)  (L-K).  (8.9) 

The  curve  in  Fig.  13  is 

r  K-(i->i;2K 

\°  =  max '  0,  - - - —  .  (8.  10) 

1  2£  (i  -n  J 

As  the  figure  labels  indicate 

|  Pjy  |  <  |  Pm  |  for  all  X<  \°  (8.  11) 

and  conversely. 

9.  ASYMPTOTIC  EFFICIENCY  USING  CARTESIAN  DATA 

The  logarithm  of  the  likelihood  of  the  n  pairs  of  observations  (1.  1),  with  a 
linear  phase  model,  is 


logL  =  -nlog7r  -  nlogfl  -  —  S 

V 


n 


t=l 


/x^  =  pcos(a>t-hp)  =  p  sin(a>t+<p)  • 


In  addition  to  (1.  3),  there  are  now  two  more  parameters: 


e4  =  w. 

We  have 


(9.  1) 


(9.  2) 
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with 


9  log  L 
9(9, 


9 


9  log  L 
99. 


S 


n 

2Z[<xt 

t=i 

9^lt 

30  +(Vt'l‘2t) 
i 

9p 

-sfi  (i=2’3’4) 

1 

(9. 3) 

Mit 

M2t 

9 

9  02 

Mlt/29A 

"  2t/2a 

9 

903 

_M2t 

"it 

9 

9  *4 

-tM2t 

•"it 

(9.  4) 

The  four  Likelinood  equations  are  obtained  by  equating  (9.  3)  to  0.  We  note  that 


n 


9.  =  ~ 
1  n 


((xt-^lt)2  +  (yt-?2t)2] 


t=l 


where  g  means  p  with  the  parameters  replaced  by  any  values  for  which  the  other  three 

equations  are  satisfied.  In  the  degenerate  case  of  constant  regression,  i.e.  w  =  0, 

u.  =  x  and  u„  =  y.  Then 
plt  2t  y 


9 


1 


is  distributed  as  — 
2n 


2 

*2(n-l) 


where,  as  always,  9  is  9^. 

Each  of  the  four  random  variables  in  (9.  3)  has  zero  expectation.  The  informa¬ 
tion  contained  in  a  sample  (x^.y^), .  . .  ,  (x^,  y^)  of  size  n  is  given  by  the  matrix 
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(9.  5) 


/  \ 

In  most  problems  we  have  J  =  n  J  ,  where  J  is  the  matrix  of  information  contained  in 

a  single  observation  having  the  form  (6.  12).  Here,  however,  we  will  find  that  the 

3 

information  contained  in  the  sample  concerning  the  frequency  oo  increases  like  n 
rather  than  n. 

The  random  variable  S  in  (9.  1)  is  distributed  as  0  •  \  ,  so  that  (see  the 


2n 


formula  preceding  (3.  5)) 


k  „  k  T  (n+k) 
(n)  - 


Thus,  after  squaring  the  first  equation  in  (9.3)  and  taking  expectations 


a1 

s  Ja=n. 


In  view  of  independence,  we  have 


It  \  2  /  3g  \  2 


90. 


+ 


2t 


dO. 


t=l 


for  i  =  2,  3  and  4.  Using  (9.  4),  and  the  formula  for  the  sum  of  the  first  n  squares, 
we  obtain 

2  (n)  =  n_ 

9  J22  2X 

92}^  =  292\n 

02J^  =  |  0“Xn(n+l)  (2n+l). 
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It  is  easy  to  see  that  the  covariance  between  the  first  and  second  line  of  (9.  3)  vanishes 
for  each  i  because  it  involves  the  first  and  third  moments  of  0  mean  normal  variables. 
That  is  to  say. 


Furthermore, 

since  the  "inner  product" 
remaining  one  we  find 

These  combine  to  give 


/n>  = 


For  the  inverse  we  have 


T(n)  _  T(n)  _  (n)  = 
J12  J13  J 14 


(n)  _  (n) 

J23  J24 


=  0 


between  the  corresponding  rows  in  (9.  4)  vanish.  For  the 


02J^  =  02  Xn(n-t-l). 


-4  0 


2  e2\ 


0  0  2Xn  Xn(n+1) 


0  0  \n(n+l)  -  X  n(n+l)  (2n+l) 


(9.6) 
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,(n)_1  _  1 


2e\ 


0 

0 

2n+l 

A(n-l) 

_ 3_ 

A(n-l) 


0 

3 


A(n-l) 

6 


A(n-l)  (n+1) 


(9.  7) 


The  marginal  asymptotic  distribution  of  9  -  9  )  and  \/ n(  0.,  -  9  )  is  normal  with 
covariance  matrix  I?  given  by  the  upper  left  hand  2x2  block  of  nj.'  .  Thus,  in 
the  notation  of  Sec.  8, 


where 


P° 


ML  (Cartesian) 


P° 


1  0 
0  2A 


(9.  8) 


This  verifies  (8.  8). 

Because  of  the  zeroes  in  (9.  6),  the  matrix  (9.  8)  is  the  covariance  matrix  of  the  ML 

estimates  of  (9.,  0j  when  0O  and  0  .  are  known.  In  fact,  for  an  arbitrary  known 
12  3  4 

phase  sequence  {(pj  in  (1.  1),  we  have  from  (9.  3)  and  (9.  4) 

n 

1  V  r/  -  .2  -  .  x2n 

9±  =  -  ^[(xt-p  cos  0t)  +  (yt  -  P  sin ] 

t=l 

n 

02  =  p2  =  [  ~  ^  (xt  cos  ^  +  yt  sin  c^.)] 2  (9.  9) 

t=l 
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10.  NUMERICAL  RESULTS  AND  RECOMMENDED  USAGE 

In  this  section  we  collect  the  results  of  desk-calculator  computations  of  the 
various  formulae  derived  in  the  preceding  sections.  On  the  basis  of  these  we  construct 
a  flow  chart  stating  which  method  is  to  be  used. 

Table  I  contains  values  of  the  asymptotic  variances  of  and  the  squared  correlation 
between 


for  Method  II  defined  by  (4.  10).  Then  entries  are  plotted  in  Figs.  2,  3,  and  4.  As 
already  mentioned,  the  corresponding  quantities  for  Method  I,  as  well  as  the  method 
of  ML,  have  the  same  appearance.  The  latter  are  relatively  difficult  to  compute,  and 
for  those  reasons  we  have  not  done  the  calculations. 

Figures  5  through  9  are  plots  of  the  entries  in  Table  2.  The  values  of  the 
generalized  variance  of  Method  I  (GVI)  were  obtained  by  means  of  the  formula 

2 

16(  1 +  A)(  1  ■+■  -4  [  ( |+  5X  +  6 A.2  +  2A3)  +  X(1  +  2X)(3  +  2X)  — 5 -  +2?C(1  +  X)] 

2  1  4  2  m 

7rm.  m .  1 

1  1  (10.  1) 

starting  with  8  significant  digits  for  ^ 

m  =  e  2  I  (^)  (v  =  0,l) 

V  V  1 

given  in  Table  9.  8  of  the  handbook  [  1  ]  .  The  arguments  X/2  for  the  *  values  are  not 
in  this  table.  These  were  computed  from  the  formula  (4.  8)  as  written,  starting  with  6 
digits  for  F,  F'  and  G  which  resulted  from  single  precision  7094  computations.  By 
comparison  with  the  minimum  4  digit  accuracy  (usually  5)  of  (10.  1),  it  appears  this 
procedure  is  satisfactory  for  plotting  purposes  up  to  X  =  10.  However,  for  X  =  16, 
only  2  good  digits  remain  and  for  X  =  20  but  1. 
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The  function  CiVI  in  Fig.  5  has  the  proper  large  argument  slope,  4A,  but  a 

graphical  intercept  of  3.  Evidently,  then  A  =  20  is  not  sufficiently  large  for  the 

asymptotic  linear  approximation  (5.6)  to  be  of  use. 

The  values  of  GVML  were  obtained  from  Benedict's  numerical  integration  of 

(B2)  to  give  1  as  a  function  of  ot  =  n/2A  . 

The  entries  in  Table  3  are  respectively  plotted  in  Figs.  10,  11  and  13. 

Table  4  gives  some  values  of  (8.4).  The  smaller  numbers  in  parentheses 

are  those  with  K  redefined  by  (8.  6).  This  reduction  in  generalized  variance  obtains 

2 

if  we  switch  from  Method  III  to  Method  IV  and  the  point  (£  ,  A)  falls  below  the  boundary 
curve  in  Fig.  13. 

The  generalized  variances  of  Methods  I  and  II  depend  on  A  =  0  /6 and  those  of 

Methods  III  and  IV  on  |  =  coscu  as  well  as  A.  If  we  knew  these  values  (which  we  don't) 

we  could  specify  which  procedure  should  be  used.  It  is  nonetheless  reasonable  to  make 

a  choice  on  the  basis  of  sample  quantities.  In  the  flow  chart  for  coherent  data  pre- 

2 

sented  below,  one  thus  uses  as  values  of  £  and  A  estimates  computed  from  the  data 

+ 

via  Method  III,  i.e.  from  (2.  9)  and  the  ratio  of  the  two  estimates  in  (2.  11)).  The 
2 

function  A  of  £  is  defined  by  (8.  10)  and  (8.  2).  The  function  K 


® 


exceeds  1  on  the  interval  (.  15, .  85).  If  such  is  the  case,  we  have  assumed  the  least 
favorable  situation,  viz.  that  K  is  maximum,  and  switched  to  Method  I  if  A  exceeds  the 

t  We  could  of  course  set  up  a  statistical  decision  procedure  using  the  large  sample 
distributions  of  the  estimates  of  ij  ^  and  A. 
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value  at  which  the  line  GVIII  x  passes  up  through  GV1  (approximately  10  as  shown 
in  Fig.  5).  From  Table  4  we  see  that  the  inclusion  of  Method  IV  gives  a  relatively 
small  improvement  over  Method  III. 

For  incoherent  data  we  use 

I  if  X  >  .  4 

II  if  X  <  .  4 

with  X  computed  as  the  ratio  0*/#*  via  (2.  4)  or  (2.  5).  (10.  3) 
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Table  1.  The  Variances  and  Squared  Correlation  for  Method  II. 

Note:  The  correlation  is  the  negative  square  root  of  the  entries 
in  the  third  column. 


X 

^11 

q22 

Q2 

M2 

^11^22 

.  1 

142.00 

141.20 

.99 

.  2 

47.00 

46.40 

.  97 

.3 

26.44 

26.04 

.94 

.4 

18.25 

18.05 

.90 

.  5 

14.00 

14.00 

.  86 

.8 

8.56 

9.  16 

.73 

1 

7.00 

8.00 

.64 

1.5 

5.  11 

7.  11 

.47 

2 

4.25 

7.  25 

.34 

2.  5 

3.76 

7.  76 

.  26 

3 

3.44 

8.44 

.  21 

4 

3.06 

10.06 

.  14 

5 

2.84 

11.84 

.  10 

6 

2.69 

13.69 

.078 

7 

2.59 

15.59 

.063 

8 

2.52 

17.52 

.052 

9 

2.46 

19.  46 

.044 

10 

2.41 

21.41 

.039 

12 

2.34 

25.34 

.033 

14 

2.29 

29.29 

.025 

16 

2.25 

33.25 

.021 

18 

2.23 

37.23 

.018 

20 

2.20 

41.  20 

.016 

24 

2.  17 

49.  17 

.013 

28 

2.  14 

57.  14 

.011 

32 

2.  13 

65.  13 

.0092 
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Table  2.  Generalized  Variances  of  Methods  1  and  II  and  Their 
Efficiencies 


A  I  II  I/II  $  ML  ML/I  ML/II 


.  1 

214.5* 

169.4 

1.  27 

.  2 

73.5 

64.8 

1.  13 

.4 

31.3 

31.9 

.981 

.5 

25.2* 

27.0 

.933 

1 

16.2 

20.0 

.810 

2 

15.3 

20.3 

.  753 

3 

17.  5 

23.  1 

.759 

4 

20.  7 

26.6 

.779 

4.  5 

22.4* 

5 

24.  2 

30.2 

.801 

6 

27.9 

34.0 

.821 

8 

35.6 

41.  8 

.853 

10 

43.  5 

49.6 

.876 

12.  5 

16 

67.3 

73.  4 

.  912 

20 

83.2 

89.3 

.932 

24 

99.2 

105.  3 

.942 

30 

123.  1 

129.2 

.952 

40 

162.9 

169.  1 

.963 

.  83836 

165.  8 

.  773 

.  979 

. 72594 

61.3 

.834 

.  946 

. 57522 

28.3 

.904 

.  887 

. 52145 

23.3 

.925 

.863 

. 35442 

15.8 

.975 

.790 

. 21316 

15.2 

.994 

.  749 

. 15101 

17.  5 

1.000 

.  758 

. 11648 

20.7 

.  778 

. 10446 

22.4 

. 094664 

24.2 

.  801 

. 079677 

27.9 

.  821 

. 060475 

35.6 

.852 

.048714 

43.5 

.877 

.039182 

53.4 

. 030753 

67.3 

.917 

.024682 

83.6 

.936 

.020613 

99.6 

.946 
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Table  3.  Functions  Associated  with  Methods  III  and  IV 


I2  A 

2  -  VaWn(|*  -|) 

K 

X° 

0 

.  500 

.500 

0 

.05 

.  744 

0 

.  1 

.257 

.  903 

0 

.  15 

1.006 

.318 

.2 

.  143 

1.071 

.  534 

.25 

1.  Ill 

.629 

.3 

.0898 

1.  133 

.674 

.375 

1.  143(max) 

.  704 

.4 

.0679 

1.  142 

.713 

.5 

.  0625(min) 

1.  125 

.  750 

.  6 

.0661 

1.095 

.823 

.  7 

.0747 

1.059 

.  974 

.8 

.0858 

1.021 

1.32 

.  85 

1.001 

1.671 

.9 

.0982 

.  982 

2.  40 

.  95 

.963 

4.611 

1 

.  Ill 

.944 

00 

t 


From  Fig. 


11,  K  =  l-^2  at  £2 


115. 


50 


Table  4.  Generalized  Variance  of  Method  III  and  Method  IV  in  Parentheses. 
Note:  2A  is  the  Generalized  Variance  of  the  Method  of  Maximum  Likelihood. 


2^  0  .1  .2  .3  .4  .5  .6  .7  .8  .9 

1  2.00  2.80  3.14  3.26  3.28  3.26  3.20  3.12  3.04  2.96 

(3.12)  (3.22)  (3.08)  (3.00)  (2.88)  (2.72)  (2.52)  (2.28) 


2 

3.50 

4.70 

5.2 

5.39 

5.  42 

3 

5.00 

6.60 

7.3 

7.52 

7.56 

4 

6.50 

8.50 

9.4 

9.65 

9.70 

5 

8.00 

10.40 

11.42 

11.  78 

11.  84 

6 

9.  50 

12.30 

13.49 

13.91 

13.98 

7 

11.00 

14.  20 

15.56 

16.04 

17.  12 

8 

12.50 

16.  10 

17.  63 

18.  17 

18.  26 

9 

14.00 

18.00 

19.  70 

20.30 

20.40 

10 

15.50 

19.  90 

21.  77 

22.43 

22.54 

20 

31.5 

39.9 

43.5 

44.  7 

49.9 

5.39 

5.30 

5.  18 

5.06 
(4.  76) 

4.94 
(4.  19) 

7.  52 

7.40 

7.  24 

7.08 

6.92 
(6.  28) 

9.65 

9.50 

9.  30 

9.  10 

8.90 
(8.  55) 

11.  78 

11.60 

11.36 

11.  12 

10.  88 

13.  91 

13.70 

13.42 

13.  14 

12.86 

16.  04 

15.80 

15.48 

15.  16 

14.  84 

18.  17 

17.  90 

17.  54 

17.  18 

16.  82 

20.30 

20.00 

19.  60 

19.20 

18.  80 

22.  43 

22.  10 

21.66 

21.  22 

20.78 

44.  7 

44.3 

42.4 

41.6 

40.  7 
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APPENDIX  A 


THE  DELTA  METHOD 


A  standard  technique  in  large  sample  distribution  theory  is  called  the  "delta 
method."  We  state  it  here  in  general  form,  since  it  is  used  throughout  the  paper. 

Suppose  t  „,t  . . .  is  a  sequence  of  estimates  of  a  p-dimensional  parameter  0 

—1  —2  — 

such  that  the  distribution  of  ^n(t  -0)  converges  as  n  —  00  to  that  of  a  p-variate 

~ n  — 

normal  random  variable  with  mean  vector  £  and  covariance  matrix  A  (not  necessarily 
of  full  rank).  We  express  this  in  symbols  by 

*s/n(t  -0)  ~  N(0,  A). 

— n  —  —  — 


Let  g.  be  a  function  of  p  variables  such  that 


3g.(x  ,  x  , . . .  ,  x  ) 

dij(-)=_J  (i  =  1,  2, . . . ,  p) 


are  continuous  at  x  =  6  with  at  least  one  nonzero  for  each  j  =  1,  2, . . . ,  q.  Put 


D  =  [dLj(£)  ’•  i  =  1>2 . P  :  j  =  l,2,  ...,q] 

and  let  g(x)  be  the  column  q-vector  with  j  component  g^(x).  Then 

N/n (£(0  -  £(6)  ~  N(0,  D'AD). 


The  conditions  of  this  theorem  place  no  requirement  on  the  existence  of  moments. 

th 

In  our  applications,  however,  the  ij  entry  of  A  is  the  limiting  value  of 

n£(t  ■ -0)  <t  -0)'. 

^  — n  —  - n  — 
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APPENDIX  B 


AN  INTEGRAL 


In  Rice's  notation, 


oi  =  p/a 


the  integral  (6,  14)  is 

4>  +  1 


I^(cbc) 

Iq(«x) 


i.  e.  ot  =  2A, 


xe  2^x  +Qf  ^(Q!X)dx. 


(B.  1) 


(B.  2) 


The  function  following  the  center  dot  is  the  density  of  x  =  r/cr,  i.  e.  a  y  '  variate.  In 
terms  of  A  we  have  the  small  and  large  argument  expansions 


1  -  2A  +  5A  +  0(A  ) 


*  =< 


1  1 
2X  '  sx2  + 


The  approximation 


as  A  —  0 

as  A  —  00  . 


(B.  3) 


1 

*  ~  1  +  2A 


(B.  4) 


agrees  with  the  result  of  numerical  integration  to  2  digits  over  the  entire  range  . 
Referring  to  (6. 15),  we  see  that  the  information  matrix  will  therefore  be  close  to 


singularity. 


t  The  author  wishes  to  thank  T.  R.  Benedict  of  the  Cornell  Aeronautical  Laboratory 
for  the  results  of  his  investigation  of  the  integral  $+1  and  the  values  obtained  by 
numerical  integration. 
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APPENDIX  C 


PROOF  OF  THE  THEOREM  OF  SECTION  7 


We  set 


n 


n 


(n)  =  1 
h  n 


h  n 


t=l 


t=l 


n 


t=l 


and  drop  the  subscript  1  from  since  for  the  present  we  are  dealing  only  with  { x^} . 
By  the  assumption  (7.  6a) 


(Cl) 


The  arithmetic  means  are  centered  at  expectations  and  the  order  term  is  a  sure  one. 
(7.  5)  is  an  almost  immediate  consequence  of  the  following  Law  of  Large  Numbers  for 
dependent  random  variables  (see  Parzen  [6]  ,  p.  419)  :  Let 

n 

1  V 

v  =  —  >  v 

n  n  t 

t=l 

where  v^,  v^,  ...  is  any  sequence  of  0  mean  random  variables  having  uniformly 
bounded  variances.  Then  as  n  •—  00 


n  n 


v  — •  0  a.  s. 
n 


if  fw  v  =0(l/n  ).  for  some  5  >  0. 


n  n 
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We  have 


n 


n-1 


f  X^F  u  |  =  |  -  N 
^  h  n  n+h  n  A 


-  )  FFa  Is  COnst' 
n  /_,  t  n  n-h  n 


a  I  =  0(l/n) 


t=l 


t=0 


by  (7.2),  and  similarly  for  Y.  For  Gaussian  processes 


f  u.u  u  ,,  u  ..  =  a.  a,  ,  +  a,  a,  ,  +  a,  a, 

L  t  t+h,  t+h. 


Thus 


and 


1  l,,*2  t+h3  hl  h2_h3  h2  hl"h3  h3  hl‘h2 


f  Z^n\u  u  -  a  )  =  -  \  ( (fu  u  u  u  -  cr^  ) 
e  h  n  n+h  h  n  ^  t  t+h  n  n+h  h 

t=l 

n 

1  Y,2 

=  —  >  (a  +  a  ,  a  ,  )  , 
n  A  n-t  n-t+h  n-t-h 


t=l 


n-1 


(C2) 


'AAw  VI  «Ft2+  IVhVh  =  0<1/n) 

t=0 

by  the  Schwarz  Inequality  and  square  summability  of  {a^} .  This  proves  the  first  statement 
in  (7.  5). 

Multiplying  (Cl)  through  by  \fn  gives 


u<n)  =  ^  x<n)  +  -Jn  Y<n>  +  Z<n>+  o(l) 

h  h  h  h 


(C3) 


Since  for  0  mean  Gaussian  processes 


r  u  u  ,,  u  ,,  =0 

c  t  t+h,  t+h. 


(C4) 


1  2 


for  all  indices,  we  have 
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n  n 


n  =  -  N  ^F(fu  u  u  • -  t  u  a  )  =  0 

L  a  b  n  / ,  / ,  t  *■  t+a  s  s+b  1  t+a  b 


t=l  s=l 


Similarly,  n  =  0  for  all  a,  b.  Thus, 


<fu (n)  uln)  =  [X(n)X(n)  +  Y^Y^l 

^  a  b  ^  ab  abJ 


+  n/[X(n>  Y<n)  +  X<n)Y<n)  ] 
^  a  b  b  a 


+  n(5z(n)Z(n)  +  o(l). 

^  a  b 


(C5) 


The  limiting  value  of  the  penultimate  is  easily  evaluated.  We  have  from  (C2) 

n  n 


„£z(n>z(n>  =  i 

^  a  b  n 


u  u  ,  u  u  -  a  o  ) 
t  t+a  s  s+b  a  b 


t=l  s=l 
n  n 


(a  a  +  a  ,,  a  )  . 

s-t  s-t+b-a  s-t+b  s-t-a 


t=l  s=l 


There  are  n-  |  k  |  indices  t,  s  =  1,  2,  . . .  ,  n  for  which  the  difference  s-t  equals 
k  (k  =  0,  ±  1,  ±2,  . . . ,  ±(n- 1)).  Hence  as  n  —  00 

Jkl 

|k|  s  n-l 


n(fzWz(n>=  ^  (‘TllYwb^Wh1 


(a  a  +  a  cr  ) 

v  k  k+b-a  k+b  k-a; 


k=-° 


°k^°k+b-a  +  °k+b+a^ 


k=-°c 


(C6) 
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r  2 

by  the  Kronecker  Lemma  and  the  summability  of  (cr  }.  (We  repeatedly  use  the  fact 

that  cr  =  <7,  to  rewrite  double  infinite  summations. ) 

-k  k 

The  main  burden  of  the  entire  derivation  of  (7.6)  is  contained  in  the  evaluation 
of  the  limit 


L(a,  b)  =  lim  n£Y(nVn) 
n  L  a  b 


n  n 


,  a 

s+b  s-t 


s=l 


(C7) 


for  arbitrary  a  and  b.  By  appropriate  changes  of  the  summation  indices  we  see  that 


lim  n(5x(n)X*n)  =  L(-a,  -b) 
n  a  b 

lim  n(5x(n)Y'(n)=  L(-a,b)  lim  n£x[n)Y(n)=  L(a,  -b)  (C8) 

nc  ab  n^ba 


since  L  is  symmetric. 

We  compute  L  by  the  following  method.  Let  m  =  m(n)  tend  to  infinity  over 
positive  integers  in  such  a  way  that 

2 

m  /n  —  0. 


Introduce 


and 


if  |  k  |  ^  m 
if  |  k  |  >  m 


n 


n 


Ln(a>  b)  =  n 


F  F  o(n) 
t+a  s+b  s-t 


t=l  s=l 


We  have 
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n  n 


|n£Y(nMn)  -  L  (a,  b)  | 
1  ^  a  b  n  ' 


const. 


n 


„  >> 
a  -a 

s-t  s-t 


t=l s=l 


=  const 


.  J  (i-iiL)|a  a( 

L,  n  1  k 


(n), 
k  1 


k  ==  n-l 


=  const. 


y 

Lj 

m+1  <  Ikl  ^n-l 


(1 - •)  |a  =  o(l) 

n  k 1 


as  n  >  m  — •  00  by  (7.  2)  and  Cauchy's  criterion.  Thus 


L(a,  b)  =  lim  L  (a,  b). 
n  n 


With  a  and  b  understood,  we  substitute  k+t  for  s  in  L  and  get 

n  ° 

n  n-t 

»(n)=  Y  f 

\  l_j  k+t+b  k 


L  =  -  \  F  ^ 
n  n  l_j  t+a  t 


(n) 


t=l 


k=-t+l 


n  will  certainly  exceed  2m  for  all  large  values  of  the  former,  and  we  then  find  by  the 

definition  of  {  of  n^} 
k 


m 


-t+1 


m 


.  (n)  . 

•T  =  < 
t 


F  a 
k+t+b  k 


k=-m 


n-t 


-m 


if  1  <  t  ^  m 


if  m+1  ^  t  ^  n-m 


if  n-m+1  ^  t  ^  n 
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Thus 


n-m 

L  =-  \  F 
n  n  L-t  t+a 

t=m+l  |  k  |  ^  m 


m 


n 


y 

/ 

Lj 


^  1 
F  cr  H — 
k-K+b  k  n 


(  'S  +  \  ^  F  T(n) 

\Lj  L  J  t+a  t 


t=l  t=n-m+l 


In  the  first  and  third  lines  in  the  preceding  display  there  are  at  most  2m  summands 

each  under  uniform  bound,  i.  e.  both  max  |  -1^  and  max  ]T^|  are  —  const,  m. 

l^t^m  n-m+l^t^n 

Thus,  after  interchanging  the  order  of  summation, 


L 

n 


n-m 

)  F  F 
/_!  t+a  t+k+b 

t=m+l 


We  split  up  the  sum  within  square  brackets  into 

n  m  n 

l-l'l 

t=l  t=l  t=n-m+l  . 


Both  the  second  and  third  of  these  is  in  absolute  value  ^  const,  m  where  "const.  "  does 
not  depend  on  k.  Hence 


L 

n 


1 0k 


kh£m 


n 


F  F 
t+a  t+k+b 


L  t=1 


=  \  a  -  \  F 

Lj  k  n  /_/  t 


k  ^m 


t+k+b -a 


t=l 


+  o(l) 


by  (7.  2)  and  the  choice  of  m.  It  is  at  this  point  that  uniformity  in  h  of  the  order  term 

hypothesized  in  (7.  3a)  is  needed.  Accordingly,  the  average  remaining  within  square 

(F) 

brackets  differs  in  absolute  value  from  A,  ,  by  a  term  o(l/  v n)  which  does  not  depend 
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on  k.  Again,  using  the  absolute  summability  of  {a  }  and  the  fact  that  m/n  —•  0,  we 
therefore  have 


as  n  —  00 .  Thus, 


L 

n 


/  ,  k  k+b-a 
|k|  —  m 


+  o(l) 


L(a,  b) 


oo 


x(F) 

k+b-a 


k=-°° 


(C9) 


for  any  fixed  integers  a  and  b. 

With  this  result,  together  with  (C6)  and  (C8),  we  return  to  (C5)  and  take  limits 
as  n  — ■  °°.  Since  X  and  X  and  a  =  a  ,  we  have  from  (C9) 

~K  K  ”  K  K 

L(-a,  -b)  =  L(a,  b)  L(a,  -b)  =  L(-a,  b). 


Consequently, 


lim  (fu 


(n)u(n)  _ 


2  1 

k=-°° 


OC 


k=-°c 


which  is  the  assertion  (7.  6a). 

The  preceding  suffices  to  prove  the  balance  of  the  theorem  with  a  minimum  of 
effort.  We  redefine  X,  Y  and  Z  by 


ii 

(n)  1  \  tj 

X'  =  -  >  F  u 
h  n  Lj  t  2, 


,(n) 


Yh  n 

t=l  t=l 


t+h  *h 

n 


li 

•11 


G  ,  u 
t+h  It 


7(n)-  iVu  u 

h  n  L  Ult  2,  t+h' 
t=l 


each  of  which  has  mean  0  because  {u^}  and  {u^}  have  0  mean  and  are  independent. 
From  (7.  3c)  we  have 
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c<n)-  HU  =  X<n>  +  Y?n)  +  Z<n>  +  o(l/N^)  . 
h  h  h  h  h 


The  conditions  of  the  Law  of  Large  Numbers  are  again  satisfied,  proving  the  second 

.  .  rv(n)  (n)_  rv(n)7(n)_  _ 

statement  in  (7.  5).  As  far  as  cross-moments  go,  not  only  is  £  -£  Y&  zfa  -  u 

but  now  also 

r(n)Y(n)  =  Q 
^  a  b 


for  all  a  and  b.  Thus, 


We  have 


f  W(n)W(n)  =  n£x(n)x£n)  +  nfY(n)y[n)+  nf  Z(n)Z<n)+  o(l). 
a  b  uab  ^  a  b  a  b 


n  n 


n^Zf>Zbn)=  zl  £<StU2,tt*Uls"2, 


s+b 


t=l  s=l 
n  n 


a  a 

s-t  s-t+b-a 


t=l  s=l 


(1 - )a,  ^  LU 

n  k  k+b-a 


k  £n-i 


akak+b-a  ' 


k=-°° 


Using  (7.  3a)  and  (7.  3b),  and  the  assumption  of  a  common  covariance  sequence,  our 
previous  evaluation  of  the  limit  L  and  (C8)  combine  to  give 
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oo 


nlim<fx<n)X<n) 
_  ^  a  b 


n  lim  ^Y(nMn)  = 

n  a  D 


k=- ' 


„  ,(F) 

k  k+b-a 


a 

k  k+b-a 


k=-8 


These  prove  (7.6b). 

We  next  have,  after  returning  the  subscript  1  in  (C3), 

(n)w(n)  =irSFu  ^  +  ZFfi  u.  +  S(u,  u  -  a  )] 

a  b  n  1  t  l,t+a  t+a  It  It  l,t+a  aJ 


l2FtU2,t+b+  SGttt>Ult  +  ZuitU2,tJ  +°(1) 


where  all  summations  are  on  t  from  1  to  n.  After  taking  expectations  we  are  left 
with  only  two  terms,  viz.  the  product  of  the  first  in  the  first  line  with  the  second  in 
the  second  plus  the  second  with  the  second.  Thus, 


n  n 


n  n 


(n)w(n)  =i  \  V  F  G  or  +  -  Y  \  p  G  a  +  o(l) 
n  Lj  Lj  t  s+b  s-t -a  n  lj  l_j  t+a  s+b  s-t 


£  u(n)wu 

^  a  b 


t=l  s=l 
n  n 


t=l  s=l 
n  n 


=  i\  Yf  G  a  +  —  "S  ^  F  G  a  +o(l). 

n  Ij  t-a  s+b  s-t  n  /_j  Lj  t+a  s+b  s-t 


t=l  s=l 


t=l  s=l 


We  return  to  (C7),  and  its  sequel,  and  everywhere  replace  Fg+b  by  Gg+b<  The  argu¬ 
ment  results  in 
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by  (7.  3c).  The  previous  two  statements  prove  (7.  6c). 

The  asymptotic  normality  follows  immediately.  Indeed,  both  the  *s/n  *s  and 

\l n  Y^’s  are  normal  for  every  n,  and  the  \ In  Z^'s  (again  for  either  definition)  are 
h  n 

known  to  be  normal  about  0  in  large  samples  with  the  covariances  of  the  limiting  dis¬ 
tribution  equal  to  the  limit  of  the  corresponding  covariances. 

Extension.  If  we  drop  the  assumption  that  {u^}  is  Gaussian,  but  retain  that  of 
summable  second  and  fourth  order  moment  sequences,  as  well  as  zero  first  and  third, 
very  little  is  changed.  The  only  place  where  (C2)  explicitly  entered  was  in  the  evalua¬ 
tion  of  lim^U^uf1^  .  In  the  non-Gaussian  case,  the  difference  between  the  two 
^  a  b 

sides  of  (C2)  is  the  fourth  order  cumulant  sequence  of  {u}. 


'*NG<hl’h2'h3)=  £Vt+h,W 

1 


2  t+h3 


1  2 


+  ah  ah 

h2  hl 


+  Vah 


The  effect  of  this  is  to  add 

oc 

1  “NG<aik’lrtb) 

k=-00 

to  the  right  side  of  (C6).  Suppose  we  write  {u^}  as  a  linear  process,  with  summable 
impulse  response  coefficients,  in  independent  identically  distributed  random  variables, 
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say  {vj, 


where 


We  can  then  evaluate  the  sum  to  get 

OO 

V  *4 

2,  *Va’k’k+b)  =  “2  ffaV 

k=-°°  K2 

2 

The  ratio  of  cumulants  k  I k~  is  the  so-called  coefficient  of  excess  of  {v  }  over  a 

T  ^  L 

N(0,  k9)  distribution,  and  it  vanishes  when  {v^}  is  so  distributed. 

In  summary,  then,  if  we  drop  the  normality  assumption,  but  retain  (C4),  then 
*4 

—  a  a  is  to  be  added  to  (7.  6a),  and  (7.  6b)  and  (7.  6c)  are  to  be  left  unaltered. 

2  a  b 
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APPENDIX  D 


THE  UNIFORM  ASYMPTOTIC  ST ATIONARITY  OF  SINUSOIDAL  REGRESSIONS 


We  may  obviously  ignore  the  multiplier  p  in  (7.  11).  Then,  using  standard 
trigonometric  additional  formulae,  for  arbitrary  integers  t  and  h 

FtF  =  [cos  cut  cos  <p  -  sin  cut  sin  (p]  [  cos  cu(t  +h)cos  cp  -  sincut(t  +h)sin  <p] 

2  2 
=  cos  <pcoscu(t+h)coscut  +  sin  <psincu(t  +h)sincut 

-  sin<p  costf[  sincu(t  +h)cos  cut  +  cos  cu(t  +  h) sin  cut] 

2  2 
=  \  cos  psf  cos  cuh+  cos(2cut  +  cuh)]  +  \  sin  <p[coscuh  -  cos(2cut  +cuh)] 

-  sin<pcos  cp  sin(2cut+  cuh)  .  (Dl) 

If  cu  is  any  multiple  of  7r,  the  second  and  third  terms  in  the  last -written  expression 
vanish  by  27r-periodicity  and  sin  mr  =  0.  The  first  gives 

n 


1 

n 


F  F  ,  =  cos  <p  cos  cuh 

t  t+h 


(cu  =0,±7r,±27r,  ...)  (D2) 


t=l 

for  every  n.  In  the  contrary  case,  we  write  (Dl)  as 

2  2 

F  F  =  icoscuh  +  i(cos  <P-sin  <P)cos(2cut  +  cuh)  -  sin</>cos  <psin(2cut+ cuh) 
t  t+h  2  2 

=  \  cos  cuh  +  geos  2<pcos(2cut  +  cuh)-  \  sin  2<psin(2cut  +  cuh) 

=  \  cos  cuh  +  |cos(2cut  +  ^ )  (D3) 


where 


ip  =  cuh  +  2  (p 
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doesn't  contain  t.  Now  we  have  the  closed  form  summations  (Knopp  [5]  ,  p.  480) 


\  .  sin(2n+l)w  ,  sin  nw  cos(n+l)cu 

i_j  2sina>  sinco 

t=l 

n 

V  „  cosco  -  cos(2n+l)ca  sin  na>  sin(n+l)w 

N  sin  2u>t  = - — -  - - r— ; - 

/_j  2sinco  sinw 

t=l 


(D4) 


Both  of  these  remain  bounded  as  n-*  ®  when  cu  is  not  a  multiple  of  n.  Thus,  after 
averaging  (D3), 


n 


-  )  F  F  =  koswh  + — — 
n  L  t  t+h 

t=l 


K 

r 

n 


(D5) 


where  K  =0(1)  depends  on  h  only  through  cos  ip  and  sin  ip  ,  and  hence  can  be 
n  h  h  /pj 

bounded  independent  of  the  lag.  This  together  with  (D2)  verifies  (7.  3a)  with  A^  given 
by  (7.  12). 

Since 

G  =  p  cos(ct>t  +  cp  -  —  ) 
t  ^ 

7T  /Q\ 

we  need  only  replace  <p  by  cp-  —  in  the  above  argument  to  obtain  the  formula  for  A^  ’  . 
For  gp  we  proceed  in  the  same  way.  Instead  of  (Dl)  we  have 

2  2 

F^G^h=4cos  cp[sin  a>h  +  sin(2o;t+aJh)]  +  h  sin  <p[  sincoh-sin(2cut  +  coh)] 

+  sin  cpcos  <pcos(2a;t  +cuh), 

and  hence 


—  ■  F  G  ,  =  sin  (ficos  cpc os  a>h  =  \  sin  2<pcos  a>h  (ca  =  0,  ±7r,  ±27r, . .  . )  (D6) 

n  l_j  t  t+h 


t=l 

for  all  n.  The  expression  corresponding  to  (D3)  is 
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F  G  =  A  sincoh  +  \  sin(2u>t  +  ij>  ) 
t  t+h  2  2  h 


with  $  unchanged.  Thus,  for  the  same  reasons, 


n 


£  2  FtGt+h  =  *  Sin  Wh  +  °(1/n)  (W  *  °’  ±7r' ±27r’  *  *  •  > 
t=l 


where  the  order  term  can  be  taken  independent  of  h. 


(D7) 
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FIGURE  CAPTIONS 


t 


Fig.  la. 
Fig.  lb. 
Fig.  2. 

Fig.  3. 

Fig.  4. 

Fig.  5. 


The  function  F  given  by  Eq.  (4.  1)  for  small  X  . 


The  function  F  given  by  Eq.  (4.  1)  for  large  X  . 

2 

The  variance  of  Method  II' s  estimate  of  2 a  given  by  Eq.  (4.  10) 
as  Qu. 

2 

The  variance  of  Method  II' s  estimate  of  p  given  by  Eq.  (4.  10) 
as  Q22. 


The  squared  correlation  between  Method  II' s  estimates  of  2 a 
and  p2  given  by  Eq.  (4.  10)  as  ^12^^11^22' 
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70 

71 

72 


73 


The  generalized  variances  of  Methods  I  and  II  given  by  Eqs.  (4.  8)  74 

and  (4. 1 1). 


Fig.  6.  The  asymptotic  efficiency  of  Method  II  relative  to  Method  I.  75 

Fig.  7.  The  integral  4>  given  by  Eq.  (6.  14).  76 

Fig.  8.  The  generalized  variance  of  the  Method  of  Maximum  Likelihood 

based  on  amplitude  data  given  by  Eq.  (6.  15). 

Fig.  9.  The  asymptotic  efficiencies  of  Method  I  and  Method  II  given  by  78 

Eqs.  (6.  20)  and  (6.  21). 

2 

Fig.  10.  X  times  the  large  sample  variance  of  the  standardized  cosine  79 

estimate  given  by  Eq.  (8.  5). 

Fig.  11.  The  function  K  given  by  Eq.  (8.2).  80 

Fig.  12.  The  largest  and  smallest  value  of  the  generalized  variance  of  81 

Method  III  given  by  Eq.  (8.  8). 

Fig.  13.  The  boundary  function  X°  given  by  Eq.  (8.  10).  82 


t  In  Figs.  2,  3,  5,  8  and  12  the  word  "variance(s)  should  be  qualified  by  after 

4  n 

division  by  4tr  . 
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Fig.  2 
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Fig.  5 
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Fig.  8 
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ADDENDUM 


ESTIMATION  OF  DOPPLER  FREQUENCY  AND  ACCELERATION 
USING  PHASE  SINES  AND  COSINES 


Our  concern  has  been  joint  estimation  of  the  component  powers  (1.  3).  Methods  I 
and  II  did  this  using  only  amplitude  data,  which  is  distributed  independently  of  the 
phase  c p  in  (1. 1).  Methods  III  and  IV  assumed  (1.  2)  to  be  the  case  and  utilized  both 
amplitudes  and  phases  (after  conversion  to  Cartesian  coordinates). 

Here  we  want  to  estimate  phase  parameters  from  the  data  p^,  p^>  . . . ,  p^.  The 
density  of  p^  is  f(p-<Pt)  for  0  —  p  —  27r  (or  any  interval  of  length  27r),  and  is  other¬ 
wise  zero.  The  function  f  is  even  about  the  origin  and  depends  only  on  p/a  >  0.  The 
expected  value  of  p  is  not  zero,  but  that  of  any  odd  27r-periodic  function  of  this 
difference  is.  As  a  result,  we  are  able  to  construct  consistent  estimates  of  the  co¬ 
efficients  in  any  polynomial  phase  model  which  use  only  the  sines  and  cosines  of  the 
phase  data.  Such  estimates  are  desirable  because  of  the  ambiguity  in  the  phase  data. 
If  this  were  resolved,  then  the  estimates  become  considerably  simpler  to  compute. 

We  consider  only  the  simplest  generalization  of  (1. 2);  viz. , 


<p  =  <p  +  cvt  +  a  — 
t  2 


(AD.  1) 


where  ca  ,  as  before,  is  the  Doppler  frequency  and  a  is  the  new  acceleration  para¬ 
meter.  The  estimate  (c a*,  o;*)  defined  in  the  next  paragraph  converges  to  (ca,  a)  with 
probability  one  as  n  — •  °°,  no  matter  what  the  unknown  value  of  p/a  >  0. 

Let  the  cosines  and  sines  of  the  independent  random  variables  p^,  p  ,  . . .  ,  p  be 
given.  Compute  the  four  averages 
n-2 

X  =  —  /  cosA^p  Y 

In  Li  pt 

t=l 

n-4 

x2  =  i  ^c°s<Apt-ti2pt+2)  Y 

t=l 


n-2 

1  V  .  .2 

L  =  n  ZSinApt 
t=l 

n-4 

r'n  I  =‘n(^Pt-^2Pt+2) 

t=l 


(AD.  2) 
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wherein  A  denotes  the  first  forward  difference  operator.  Graphically,  the  estimates 
are  obtained  from  D 


where  P.  =  (X.,  Y  ).  If  necessary,  360°  is  to  be  added  to  the  difference  to  put  w*  be¬ 
tween  0°  and  360°.  The  limiting  amplitudes  are  rather  complicated  homogeneous 

2  2 

functions  of  the  noncentral  parameter  X  =  p  /la  . 

If  ot  is  known,  we  estimate  a;  using 

n-1  n-1 

X  =  -  cos[Ap  -a(t+5>]  Y  =  —  /  sin[Ap  -a(t+i)]  .  (AD.  3) 

0  n  /  y  t  u  W  l_j  l 

t=l  t=l 

The  point  (X^,  Y^)  converges  to  (A  cos  w,  Asinw)  where  A  =  A(A). 

The  computation  of  PQ  and  ?1  in  terms  of  the  cospt's  and  sinp^s  is  a  simple 
task  relative  to  computing  P  .  For  the  last,  we  must  use  the  fact  that 

cos  kp  =  T^cos  p) 

sin  kp=  sinp  ^(cosp), 

where  T  and  U  are  the  first  and  second  kind  Tchebichev  polynomials  of  degree  k. 
k  k 
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A  certain  amount  of  recursive  calculation  is  possible. 

Phase  density  and  trigonometric  moments.  We  first  derive  the  density  of  with 

ip  left  arbitrary.  As  in  Sec.  3  we  take  the  bivariate  circular  Gaussian  distribution, 
t  2 

centered  on  the  point  (p  cos  </^,  p  sin</^)  with  variance  cr  ,  and  make  the  change  ol 

variables  to  polar  coordinates.  After  appropriately  completing  the  square,  the  proba¬ 
bility  density  function  of  p^  becomes 


2  2 

-2  a  sin  (p-<p  ) 


£ 


-2  [Z- 


a  cos(p-<?t)] 


ze 


dz 


Here  we  have  used 


a  =  p/a 


for  the  Rice  parameter  to  avoid  confusion  with  the  acceleration  (ci.  B.  1).  Alter  adding 
in  and  subtracting  out  a  cos(p-<?t)  from  the  integrand,  and  then  integrating  by  parts, 
there  results 

ft(p)=  f(P-<Pt)  (0  —  p  —  2tt) 


m  = 


2n 


i  2 

-2a 

e  + 


a  cos  0 e 


i  2  .  2. 

-2 a  sin  0 


a  cos  0 


1 


a,2 

2  t 


dt 


(AD.  4) 


When  a,  i.e.p,  vanishes  this  properly  reduces  to  a  uniform  distribution  over  (0,  27r). 
From  (AD.  4)  it  is  clear  that  the  random  variables 


u^  =  cos(pt-^t)  =  sin(pt-^t)  (AD.  5) 

have  the  same  joint  distribution  for  all  t.  We  will  need  their  first  and  second  moments 
in  order  to  compute  the  limiting  values  of  die  estimates. 

Letting  u  and  v  denote  the  common  variates  we  have 
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r2*-* t 

-<Pt) f.(p)dp  =  J  cos  6  f(0)d0  = 
'*t 


>27 r 

I  cos  0  f(0)d0 
0 


because  the  integrand  is  periodic  in  9  with  period  2tt  .  Since  cos  6  times  the  leading 
term  in  (AD.  4)  integrates  to  zero, 


wherein 


2  2 
a  sin  0 


G(a  cos  6)d0 


G(x)  = 


dt 


is  V2tt  times  the  standardized  normal  distribution  function.  We  split  (fu  into  the 
integral  from  0  to  7t  plus  7t  to  27r,  and  replace  6  by  6+n  in  the  latter.  This  changes 
the  sign  of  cos  6  and  sin  6  : 


But 


6e 


1 

2 


2  .  2 
a  sin 


9 


G(  -a  cos  0)d0. 


G(-x)  =  \l2n  -  G(x) 


for  all  x,  so 


<f»- 

^  27T 


n/  27T 


e 


de 


e 


d6 
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where 


A  =  a  /2 

2 

is  the  noncentral  parameter.  Substituting  sin  0  =  i(l-cos20)  and  then 
(see  (3.  2)) 

f'-Xsin2*  1  ’I  f2*^08*  „ 

\  e  =  -  e  \  e  d<p  =  ire  I  (-) 

J0  2  J0 


Since  the  derivative  of  I  is  1^,  the  preceding  combine  to  give 


r  _  n/ttA  2  .  .X.  ,  A 

u  e  [  IPi(  O)“i"l/|(o)  J 


O' 2 7  1  2 


for  the  mean  value  of  cos(p  -<p  ). 

t  2 

We  compute  the  mean  value  of  sin  (p^-^)  by  the  following  device. 


=  1. 


2tt  _i  2  2 

sin  0  cosO  e  ~a  Sin  G(a  cos0)d0 


0 

2tr 


K 


i  2  2 

/i  ”  2  s  sin.  0  .  _  v  ,  q 

cos  6  e  G(acos0)d0  , 


0 


which  we  view  as  functions  of  a.  Then 


r*2TT  o  .  _i  2 

(fv2  =  J  sin20f(0)de=  |e'2a  +  |pj 


p27t 

=  J  f( 

J0 


1  2  2 

2  a  .a 


1  =  \  f(0)d0  =  e  +  ^  K  . 


If  we  differentiate  the  integral  representation  for  K  we  get 


9  =  <p/2  gives 


(AD.  6) 


Define 
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dK 

da 


-aj 


=  27T 


The  derivative  of  the  other  expression  is 


Thus 


(AD.  7) 


after  equating  the  two  formulae. 

Any  odd  power  of  v  is  orthogonal  to  any  function  of  u,  and  hence  they  are  un- 
correlated.  We  have,  for  k=0,l,2,..., 


c  2k+l 

Cv 


1  f 


27T 


g(u)  =  2  J  sin  0[  sin 
-27T 


2k 


6  g(cos  0)f(0)]d# 


because  the  integrand  has  period  2ir.  Since  the  function  within  square  brackets  is 
even  the  integrand  is  an  odd  function,  so 


(fv2k+1g(u)  =  0. 

In  particular, 

(fv  =  0  (fuv  =  0. 


(AD.  8) 


(AD.  9) 


(u,  v)  is  an  example  of  a  pair  of  uncorrelated,  yet  obviously  dependent,  random 
variables. 
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Estimation  of  frequency  when  acceleration  is  known.  The  method  lor  analyzing 
the  large  sample  behavior  of  the  points  P  (i  =  0, 1,  2)  is  the  same  in  each  case,  but 
gets  more  complicated  with  increasing  i.  We  will  therefore  start  with  the  averages 
(AD.  3).  We  have 


Apt  =  ca  +  a(t+-2‘)  +  A(pt-</>t) 

2  2 
A  pt=  a  +  A  (pt-  <jPt)  , 

only  the  first  of  these  being  needed  at  the  present.  Setting 

n-1 

xo  4  I  cosA(>VV  Yo  4  1  sinA(  W  • 

t=l  t=l 


we  have 


XQ  =  X^cosca  -  Y^sin  ca 
Yq  =  X^sinca  +  Y^cos  ca  . 


In  the  notation  (AD.  5) 


(AD.  10) 


(AD.  11) 


cosA(pt-?t)=  Vt+1  +  Vt+1 

(AD.  12) 

sii,A(pt-«>t)  -  utv(+1  -  v(ut+1  . 

Since  (u  ,  v  )  are  independent  and  identically  distributed  we  have,  using  (AD.  6)  and 
(AD.  9) 

(fcosA(pt-?t)  =  A 
(fsinA(pt-<^t)  =  0  . 


i 
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We  have  set 


A  =  (fu  = 
B  =  (fu2  = 


7rA  -A 

Te 


vl> + 1£> 


l- 


1  -e 


-A 


2A 


(AD.  13) 


the  latter  will  be  used  subsequently.  If  we  subtract  the  constant  A  from  we  ob¬ 
tain  the  arithmetic  mean  of  first-order  dependent,  zero  mean,  bounded  random 
variables.  The  same  is  true  of  as  it  stands.  A  trivial  application  of  the  Law  of 
Large  Numbers  for  dependent  random  variables  (used  in  Appendix  C)  shows  that 
(Xq,  Y^)  converges  to  the  point  (A,0)  as  n  —  00  with  probability  one  and  in  mean 
square.  Thus,  in  both  these  senses, 


Xq  — •  A  cos  oo  Yq  —  A  sin  oo  .  (AD.  14) 

As  claimed,  the  phase  of  the  point  consistently  estimates  the  unknown  oo.  The 
function  A(A)  increases  steadily  from  A(0)  =  0  to  A(+»)=  1.  As  we  would  expect, 
estimation  is  difficult  when  A  is  small. 

The  large  sample  distribution  of  P^.  For  the  second  moments  of  the  random 
variables  in  (AD.  12)  we  find,  in  the  notation  (AD.  13), 


2  2 

B  +  (1-B) 

for 

s=  t 

(fcos  A(p  -Op  )cos  A(p  -(p  )  =  < 

L  L  u  o 

AB 

for 

|s-t  | 

=  1 

A2 

\ 

for 

|s-t| 

2=2 

2B(1-B) 

for 

s  =  t 

(fsinA(p  -<p  )sinA(p  -<p  )  =  < 

L  L  O  u 

-A(l-B) 

for 

|s-t| 

=  1 

0 

for 

|s-t| 

2=  2 

r  cosA(p  -cp  )sinA(p  -op  )  =  0  for  all  t,  s. 
^  t  t  s  s 
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We  have,  after  replacing  n-1  by  n  in  the  upper  limit  of  the  sum  which  is  , 

n  n 

n(f(Xj)-  A)2=  ^  ^  ^  (f[cosA(pt-<Pt)-Al  [ cos  A(pg-pg)-A] 


t=l  s=l 
n  n 


t=l  s=l 


We  have  set  a  =  (fcos  A(p  -<p  )cos  A(p  -<p  J  -  cr  There  are  n-  |k  |  indices 
t  -S  t  t  S  b  b-l 

t,  s  =  1,2 . n  for  which  t-s  =  k  (k  =  0,  ±1, . . . ,  =Kn-l)).  Hence 


n(f(X;-A)2=  ^  (‘•V-)(‘Tk-A2) 

I k |  -  n-1 

=  (a0-A2)+  2(1-  £)  (a4-  A2)  , 

and 

lim  n(f'(X/  -  A)2  =  B2  +  (1-B)  -  A  +  2(AB-A  )  . 

n  u 

Similarly,  we  find 

limn(fY'2=  2(1 -B)  (B-A) 

n 


Furthermore 

<f<xo'A)Yo=0 


for  every  n. 

By  Liapounov’s Central  Limit  Theorem  it  follows,  for  n  —  °°,  that  ^n(X^  -  A) 
and  's/n  tend  to  bivariate  normality  about  (0,  0)with  covariance  matrix  given  by 
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the  above  limiting  values  of  the  corresponding  second  moments.  According  to  (AD.  11), 
Vn(X^-Acos  w)  and  \Zn(Yg-  A  sinca)  are  obtained  from  these  variates  by  a  rotation 
through  the  angle  c a.  Thus,  introducing 

C  =  (1-B)2  +  (B-A)  (B+3A) 

C2  =  2(B-A)  (1-B), 

we  see  tliat  \ln{X  -  A  cos  ca)  and  'Vn^- A  sinca)  have  a  large  sample  normal  distribu¬ 
tion  about  the  origin  with  covariance  matrix 


2  2 

C  cos  ca  +  C  sin  ca  (C  -C  )sin  ca  cos  ca 

1  Z  x.  A 


2  2 

(C  -C  ) sinca  cos  ca  C  sin  ca  +  C  cos  (a 

1  Z  1  z 


(AD.  15) 


From  (AD.  13)  we  have  0  <  A  <  B  <  1  ,  so  the  C's  are  positive  functions  of  X. 

Joint  estimation  of  frequency  and  acceleration.  We  now  consider  the  point  in 
(AD.  2).  If  we  set 

n-2  n-2 

xi =  Yi =  ^  I 

t=l  t=l 


then,  according  to  (AD.  10),  (AD.  11)  is  true  with  ca  replaced  by  a  and  the  0  sub¬ 
scripts  replaced  by  l's.  The  expression  comparable  with  (AD- 12)  is 

cosd?(p  -*  )  =  ut(2uf+1-l)ut+2  +  2vtut+1vt+1ut+2 


'Vt(2uI+l',|vt+2  +  2ut"ttl'ttlV2 

sin  a2(Pt  -(0t)  =  vt(2u2+1  -l)ut+2  -  2utut+1vt+1vt+2 


(AD.  16) 


+  ut(2ut+r1)vt+i +  2,,ViViV2  • 
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From  these  and  (AD.  6)  -  (AD. 8  )  we  find 


2 

(fcosA2(pt-(?t)  =  t  u(2(fu2-l)  =  A(2B-1) 
(f  sin  A  (pt-<^t)  =  0 

Thus,  for  the  same  reasons  and  in  the  same  senses, 


(AD.  17) 


X  —  A(2B-l)cos  oi  A(2B-l)sin  a  .  (AD.  18) 

The  function  B(A)  in  (AD.  13)  increases  steadily  from  B(0)  =  1/2  to  B(+  °°)  =  1,  so  the 
limiting  amplitude  is  positive  for  all  positive 

Turning  now  to  the  second  point  in  (AD.  2),  we  have  from  (AD.  10) 


2  2 

Apt  -  t  A  pt  =  w  +  |a  +  A(pt  -<Pt)  -  t  A  (pt+2  -  <Pt+2) 

=  ca  +  5  a  +  z  (say). 

Again,  (AD.  11)  holds  with  ca  replaced  by  ca+  2  a  and  the  0  subscripts  replaced  by 
2's  when  we  set 

n-4  n-4 

X'  =  —  /  cos  z  Y'  =  —  /  sin  z  . 

2  n  A  t  2  n  Lj  t 

t=l  t=l 

2 

Since  A(p  -  )  and  A  (p  -  are  independent  random  variables, 

f  cos  z  =  A  (fcostA2(p  -  <p  ,„) 

^  t  u  Vh  t+2  t+2 

2 

(fsinzt  =  A  (f  sint  A  (Pt+2‘  ?t+9) 

results  from  the  display  preceding  (AD.  13).  Let  us  now  make  the  following  redefini¬ 
tions  in  the  right  side  of  (AD.  16)  : 


94 


"t  =  cos  t(pt+2'  ’W 


u,+l  =  cos  t(pt+3-  <Pt+3) 


V2=COSt(Pt+4-V4)  ' 


(AD.  19) 


and  similarly  for  the  v's  with  sines  in  place  of  the  cosines.  From  the  second  display 

2 

in  (AD.  16)  it  is  clear  that  the  expectation  of  sintA  (p^-  vanish  if  that  ol 

sin  t(p  -  <p  )  does  for  all  j.  Such  is  the  case: 
t+j  ^t+r 


•n+j 


p  J  p  27T 

(5 sin  t(p  (p  )=\  Sint0f(0)d0=  \  sint<?f(0)d0  =  0 
t+J  t+J  J  „  J  Q 


because  t  is  an  integer  and  the  periodic  integrand  is  odd.  (If  we  had  averaged  the 

2 

cosine  and  sine  of  Ap^-(t  +  g)A  p^  in  an  attempt  to  get  w  directly,  we  would  have  had 
a  t+g  in  the  integrand,  making  it  a -periodic.  (5 sinz^  would  therefore  not  be  zero. ) 
In  the  same  way  we  have  from  the  first  display  in  (AD.  16) 

<fcos  ta2(Pt+2-  <Pt+2>  =  (fut(2<fu2+1-l)  <fut+2  =  At(2Bt-1) 


where  the  u's  are  defined  in  (AD.  19).  (fu  and  (fut+2  are  ecluaf>  anh  we  have  set 


l  p2' 

i2  =  \ 

1  J0 


277 

cos  t0  f(0)d0 


27 r 


B.  =  ^  cos  t0  f(0)d0  . 
0 


(AD.  20) 


When  t  =  1,  these  become  the  A  and  B  of  (AD.  13). 

The  integrals  (AD.  20)  are  hard  to  evaluate.  However,  only  the  first  need  be 
considered  because 


2v4  ■  4 
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lor  all  t.  (Aj_  is  the  value  of  die  characteristic  function  at  the  point  k 
variable  p  -  )  if  we  assume  the  sequence  of  the  A's  is  such  that  the 


lim  - 
n  " 


1 

A  A;  =  L 
t  2t 


t=l 


exists  (which  is  not  an  unreasonable  supposition),  then  it  is  true  that 


’  AL  cos(cc+ l  o) 


Y7  —  AL  sin(cu+  2  a) 


asn-°° 


of  the  random 
time  average 

(AD.  21) 


(AD.  22) 
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